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Abstract 

Results  are  presented  of  a  study  to  determine  appropriate  methods  of  flow 
modification  for  mifigoting  ice  jam  formation  in  navigoble  rivers.  Based  on  a 
review  of  alternative  methods  for  flow  modification,  it  is  concluded  that  for  many 
rivers,  especially  large  ones,  the  most  appropriate  method  involves  controlled 
ice-cover  formation  through  the  regulation  of  river  flow.  Flow  discharge  and 
stoges  would  be  regulated  by  controlling  the  flow  releases  from  reservoirs  and 
flow  stoges  at  river  dams,  such  that  optimal  flow  conditions  prevoil  for  rapid 
formation,  and  subsequent  maintenance,  of  an  accumulation  ice  cover  over 
river  reaches  in  which  potentially  large  amounts  of  frozil  ice  may  grow. 
Accumulation  covers  would  be  formed  of  frazil  ice  pons  and  floes  and,  if 
appropriate,  broken  ice  conveyed  from  upstream.  Existing  dams,  augmented 
where  needed  by  navigable  ice  booms,  could  serve  as  retention  structures  for 
the  development  of  accumulation  covers.  A  preliminary  indication  of  the 
feasibility  of  this  method  for  controlling  ice-cover  formation  on  stage-reguloted 
pools  of  the  Ohio  River  is  assessed  through  the  use  of  a  numerical  model  that 
simulates  ice-cover  formation  from  frazil  ice.  It  is  found  that  this  approach  holds 
promise  for  mitigating  jam  occurrence,  although  its  implementation  necessarily 
entails  monogement  of  flow  through  major  portions  of  the  Ohio  River.  The  results 
of  fhe  study  are,  to  a  limited  extent,  generalized  to  other  rivers. 


Cover:  Severe  ice  accumulation  upstream  of  Locks  and  Dam  20,  Mississippi 
River. 
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NOMENCLATURE 


A  flow  area 

B  top  '.vidth  of  water  surface 

C  frazil  ice  concentration 

Cp  specific  heat  of  water 

d  flow  depth 

e  porosity  of  individual  ice  floes 

€(.  overall  porosity  of  ice  accumulation 

=  ep  +  (l  -ep)e 

Cp  porosity  of  ice  accumulation  between  ice 
floes 

longitudinal  diffusivity  coefficient 
F  Froude  number 

Fmax  maximum  Froude  number  for  ice-cover 
progression 

F(.  critical  Froude  number  for  juxtaposition 

g  gravitational  acceleration 

H  z  +  d+r\ 

Aja  energy  exchange  coefficient  at  air/ice 

interface 

energy  exchange  coefficient  at  water/air 
interface 

Awi  energy  exchange  coefficient  at  water/ice 

interface 

k  thermal  conductivity  of  ice 

L  latent  heat  of  fusion  of  ice 

i  channel  length 

rtf.  composite  Manning’s  coefficient 

nj,  Manning’ s  coefficient  for  channel  bed 

fij  Manning’ s  coefficient  for  ice  cover 

Q  river  discharge 

Qo  optimum  discharge  associated  with  minimum 

time  for  ice-cover  formation 

Cn  regulated  discharge  at  upstream  end  of  pool 
R  hydraulic  radius 


T  water  temperature 

Fa  air  temperature 

Fjo  constant  air  temperature 

Tf  freezing  point  of  water 

F„  temperature  of  inflow  water  at  upstream  end 
of  pool 

Fq  initial  water  temperature 

Fj  temperature  at  the  surface  of  the  ice  cover 
t  time 

/a  lime  for  an  ice  cover  to  form  over  a  pool 
Sf  friction  slope 

Sj  specific  gravity  of  ice 

V  flow  velocity 

V  volume  of  ice 

X  distance 

Yf  regulated  depth  of  flow  at  downstream  end 

of  pool 

y  z  +  d 

z  bed  elevation 

« 

jj  spatial  weighting  factor  for  the  Preissmann’s 

scheme 

E  Awa^CPjFjd) 

T\  equivalent  thickness  of  ice  cover  =  Sifl 

ice-cover  thickness 

A,,  increment  of  ice-cover  thickness 

T|j  thickness  of  ice  floe 

6  temporal  weighting  factor  for  the 

Preissmann’s  scheme 

p  density  of  water 

Pi  density  of  ice 

(!>vya  net  heat  flux  from  water  to  air,  for  unit 
surface  area  of  flow 
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Flow  Regulation  For  Controlled  River-Ice  Formation 

SUBASH  C.  JAIN,  ROBERT  ETTEMA.  AND  INBO  PARK 


INTRODUCTION 

Although  many  northern  rivers  are  regulated  for  the 
purposes  of  flood  control,  navigation,  and  hydropower 
generation,  few  are  regulated  to  mitigate  problems 
arising  from  ice  formation.  In  some  respects  this  is  not 
surprising,  as  flow  regulation  to  mitigate  ice  problems 
entails  an  additional  tier  of  constraints  on  those  imposed 
for  ice-ffee  flows.  Nevertheless,  for  rivers  that  experi¬ 
ence  frigid  winters,  the  penalties  of  not  regulating  to 
mitigate  ice  problems  can  be  severe.  The  upper  portion 
of  the  Ohio  River  is  one  such  river.  It  is  prone  to  severe 
ice  jams  that  are  of  potentially  major  economic  conse¬ 
quence. 

This  study  examines  the  technical  (hydraulic)  feasi¬ 
bility  of  regulating  river  flow  to  control  ice-cover  for¬ 
mation  so  that  ice  jams  do  not  occur.  It  was  conducted 
within  the  context  of  the  U.S.  Army  Corps  of  Engineers 
(COE)  River  IceManagement  (RIM)  Program,  which  is 
aimed  at  developingengineeringmethodsformanaging 
navigable  rivers  that  become  ice-covered  during  winter. 
Although  it  is  primarily  focused  on  ice  formation  over 
stage-regulated  pools  in  the  upper  portion  of  the  Ohio 
River,  the  results  from  this  study  are  generalized  to 
otlier  rivers. 

Background 

Ice-cover  formation  on  rivers  characteristically  in¬ 
volves  two  concurrent  processes.  One  is  the  genesis  of 
frazil  ice  from  individual  crystals  formed  in  super¬ 
cooled  water  to  ice  floes  and  ultimately  to  a  solid  ice 
cover.  The  other  is  static-  or  border-ice  growth  in  zones 
of  relatively  quiescent  flow,  notably  along  banks  and  in 
backwater  areas.  Frazil  ice  typically  accounts  for  the 
greater  part  of  initial  ice-cover  volume,  but  border  ice 
growth  plays  an  important  role  in  modifying  flow 
surface  geometry  by  narrowing  it  in  such  a  way  that,  at 
some  critically  congested  section,  frazil  floes  become 
lodged  and  in. date  upstream  progression  of  a  full,  solid 
ice  cover. 


Often,  ice-cover  formation  on  rivers  is  remarkably 
orderly,  causing  little  disruption  to  the  water  flow. 
River  water  under  frigid  air  becomes  supercooled  to  a 
few  hundredths  of  a  degree  centigrade  (or  Fahrenheit) 
and  is  seeded  by  ice  fragments  from  several  sources  so 
that  small  ice  crystals  appear.  These  crystals,  called 
frazil  ice,  grow  as  they  are  conveyed  by  the  water  flow. 
After  undergoing  several  morphological  changes,  from 
floes  to  slush,  which  rises  to  produce  ice  pans  that  float 
at  the  water  surface,  frazil  ultimately  comes  to  rest  as 
floes  of  fused  pans  lodged  in  an  ice  cover  formed  as  a 
more  or  less  single  layer  of  floes  accumulated  over  a 
river  reach.  Typically,  a  cover  begins  at,  and  progresses 
upstream  from,  a  section  critically  narrowed  by  border 
ice  growth.  Figure  1  illustrates  ice-coverformationover 
the  Cedar  River  in  Iowa.  The  physical  properties  and 
growth  cycle  of  frazil  ice  are  described  comprehensive¬ 
ly  in  several  publications,  including  those  by  Ashton 
(1986),  Ettema  et  al.  (1984),  Daly  (1984a),  Martin 
(1981),  and  Osterkamp  ( 1 978). 

Certain  flow  conditions  may  preclude  the  orderly 
development  of  an  ice  cover  described  above .  Instead  of 
coming  to  rest  as  an  accumulation  cover  of  floes  and 
pans  juxtaposed  in  a  single  layer,  fr  izil  slush  and  pans 
transported  at  relatively  high  velocities  may  be  swept 
downstream  to  form  a  jam  at  some  location  where  the 
channel  capacity  to  convey  ice  is  overwhelmed.  The  re¬ 
sult  is  a  type  of  freeze-up  jam  (lAHR  1986).  It  may 
develop  as  a  so-called  hanging  dam  beneath  an  existing 
ice  cover  if  the  frazil  is  predominantly  slush.  Alternate¬ 
ly,  if  the  frazil  is  in  the  form  of  pans,  slush,  some  floes, 
and  broken  border  ice,  a  jam  may  develop  that  thickens 
and  compacts  by  shoving  and  collapse.  Such  a  jam  is 
shown  in  Figure  2.  Whichever  type  of  jam  occurs,  chan¬ 
nels  become  impassable  to  river  traffic  and  river  flow 
becomes  restricted,  causing  a  significant  rise  in  stage 
that  usually  results  in  flooding. 

The  magnitude  and  steadiness  of  water  flow  have  an 
important  bearing  on  whether  an  ice  cover  or  an  ice  jam 
forms.  The  influence  on  ice-cover  formation  of  flow 
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Figure  2.  Frazil-generated  ice  jam  impeding  navigation  in  the  Illinois  Waterway. 


magnitude  is  customarily  described  in  ternis  of  cither 
surface  velocity  of  flow,  or  a  stability  criterion  based 
on  flow  Froude  number,  F  =  VKgdP-^.  in  which  I  and  d 
are  the  mean  velocity  and  depth  of  flow,  respectively, 
and  g  is  the  gravity  acceleration.  It  is  commonly  held 
(e.g..  Michel  1978.  Ashton  1986)  that  if  V  ^  is  less  than 
about  0.6  m/s  for  rivers  of  medium  and  large  depth,  flow 
i,s  insufficiently  agitated  so  that,  in  calm  frigid  air,  ice 
crystals  grow'ing  on  the  water  surface  directly  fomi  a 
solid  ice  cover;  for  higher  values  of  V,  ice-cover  fonna- 
tion  is  as  outlined  above.  Field  and  laboratory  observa¬ 
tions  (Kivisild  1959,  Uzuner  and  Kennedy  1972)  indi¬ 
cate  that  when  F  is  less  than  about  0.08  to  0. 1 5  (the  larger 
value  corresponding  to  deep  flows;  d  in  excess  of  about 
8  m).  individual  ice  floes  and  frazil  pans  gather  at  a 
downstream  barrieras  a  more  or  less  si  ngle  layer  that  de¬ 
velops  upstream,  forming  a  so-called  accumulation  cov¬ 
er  of  juxtaposed  pans  or  floes,  as  described  above,  Shen 
et  al.  (1984)  suggest  that,  for  this  condition,  F  =  0.05. 
Thickened  ice  covers  and  ice  jams  are  attributed  to  flows 
with  higher  values  of  F.  Large  discharges  may  preclude 
ice-cover  formation  and  contribute  to  ice-cover  breakup. 
Fluctuations  in  flow  stage  may  destabilize  ice  covers  or 
inhibit  their  formation.  An  abrupt  rise  in  stage  may 
create  uplift  forces  that  lead  to  ice  sheet  breaking  up. 

The  calamitous  consequences  of  jam  fonnation  w'crc 
dramatically  apparent  in  January  1978  at  Markland 
Locks  and  Dam  on  the  Ohio  River.  The  “Markland  inci¬ 
dent,”  as  it  has  come  to  be  called,  resulted  from  a  se¬ 
quence  of  events  in  which  a  period  of  frigid  weather, 
resulting  in  considerable  ice  production,  was  followed 
by  a  period  of  rain  that  led  to  a  fivefold  increase  in  dis¬ 


charge  over  4  days.  Ice  from  frazil  growth  and  the 
break  up  of  border  ice  was  coti  vexed  and  accumulated  as 
a  large  jam  in  the  pool  immediately  u,)sire;im  of  Mark- 
land  Locks  and  Dam.  Increasing  water  discharge  and 
stage  caused  the  jam  to  break  and  be  swept  dow  nsirc  im 
against  the  dam.  w  here  it  relonned.  causing  much  dam¬ 
age,  destroying  portions  of  the  dam  and  sinking  se\  eral 
vessels.  The  Markland  incident  and  the  ice  conditions 
le;iding  to  it  are  described  in  a  detailed  report  prepared 
by  the  Ohio  Riser  Division  Ice  Committee  (ORDIC 
1978).  It  raised  mans  issues  concerning  appropritile 
operational  and  structural  solutions  to  managing  river 
ice  prohlemson  regulated  ris  ers.  In  a  concerted  eU'  '  :  to 
re.solve  these  issues,  the  U.S.  Anns  Corps  of  Lngineers 
initiated  the  River  Ice  Management  (RIM)  Program. 

Scope  of  the  study 

The  general  objective  of  this  study  was  to  identify 
means  of  hydraulically  modifying  How  m  the  Ohio 
River  so  as  to  control  frazil  ice  formation  and  thereby 
prevent  thcoccurrenceof  ice  jams.  The  nature  of  frazil 
ice  formation  and  the  .scale  at  svhich  it  occurs  in  large, 
navigable  ris'crs  such  as  the  Ohio,  however,  makes  it 
readily  apparent  that  flow  modification  necessarily 
implies  regulatioti  of  llosv  discharge  and  stage  in  (rrder 
to  create  flosv  conditiotis  that  facilitate  accumulation- 
cover  formation  overcxletisive  lengths  of  channel  in  as 
speedy  a  tnantier  as  is  practically  possible.  In  this  cmi- 
neclion,  it  is  also  apparent  that  structural  means  ol  llosv 
modification  used  without  flosv  regulation  svould  not. 
under  ttiost  circumstances,  be  adequate  to  control  fia/il 
ice  formation. 
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Consequently,  this  study  is  specifically  aimed  at  ob¬ 
taining  preliminary  indications  as  to  the  feasibility  of 
regulati  ng  flow  for  the  purpose  of  facil  itati  ng  control  led 
ice-cover  formation  on  stage-regulated  pools  of  the 
Ohio  River  and  rivers  generally.  The  feasibility  of  con¬ 
trolled  ice-cover  formation  is  assessed  through  the  use 
of  a  computational  model  that  simulates  one-dimen¬ 
sional  flow,  frazil  ice  generation,  and  accumulation- 
cover  formation  on  pools  of  known  geometry  and  slope. 
It  is  assumed  that  accumulation  covers  form  and  progress 
upstream  from  existing  dams  or  ice  booms  that  act  as 
retention  structures.  Feasibility  assessment  is  based  on 
the  following  criteria,  which  have  to  be  considered  at 
different  locations,  or  regulated  reaches,  of  a  river; 

•  Whether  a  required  discharge  can  be  maintained 
for  the  time  needed  to  form  an  ice  cover, 

•  Whether  flow  stage  restrictions  can  be  adhered  to, 
and 

•  Whether  flow  regulation  creates  any  other  prob¬ 
lems,  such  as  channel  scour. 

An  assessment  of  feasibility  is  made  for  controlled 
ice-cover  formation  on  two  stage-regulated  pools  of  the 
Ohio  River — the  pools  upstream  of  the  lock  and  dam 
installations  at  Hannibal  and  Montgomery.  As  the  flow 
stages  for  these  pools  are  narrowly  constrained  to  within 
1 .5  m  (4  to  5  ft),  this  assessment  concentrates  primarily 
on  discharge  regulation. 

The  numerical  model  yields  information  on  periods 
of  time  for  accumulation  covers  to  form,  volumes  of  ice 
grown,  ice-cover  profile,  and  water  stage  profiles.  In¬ 
formation  on  ice-cover  formation  times  and  ice  vol¬ 
umes  is  presented  in  a  generalized  format  so  that  it  can 
be  applied  in  assessing  the  feasibility  of  controlled  ice- 
cover  formation  for  a  range  of  rivers. 

Flow  modification  to  prevent  ice  jams  that  result 
from  ice-cover  breakup  (break-upjams)  (I  AHR  1 986)  is 
not  examined  here,  except  insofar  as  it  affects  the  main¬ 
tenance  of  accumulation  covers  once  they  are  formed. 
In  some  respects,  flow  regulation  to  prevent  ice-cover 
breakup  involves  fe werconstraints  than  are  required  for 
ice-cover  formation.  Essentially,  it  entails  attenuation 
of  runoff  discharges  and  the  prevention  of  sudden  stage 
rises  so  that  ice  covers  do  not  begin  breaking  at  their 
upstream  end,  which  leads  to  a  situation  v/here  ice  dis¬ 
charge  increases  as  the  ice-cover  breakup  progresses 
downstream.  If  an  ice  cover  is  to  break  up,  it  is  best  that 
breaking  commence  at  the  downstream  end  so  that  the 
ice  discharge  does  not  cumulatively  increa.se  and  re¬ 
mains  within  the  river’s  capacity  to  convey  it.  Where 
needed,  an  ice-breaking  vessel  can  be  used  to  instigate 
ice-cover  breaking  from  downstream.  This  practice  is 
described  by  Aleinikov  and  Korenkov  ( 1 972)  and  Rozs- 
nyoi  (1972)  for  ice-cover  breakup  on  Soviet  and  Hun¬ 
garian  rivers  and  by  Petkovic  et  al.  ( 1 984)  for  the  Dan¬ 


ube  River  in  Yugoslavia.  Deck  (1984)  cites  the  use  of 
amphibious  craft  to  break  ice  on  Cazenovia  Creek, 
which  flows  into  Lake  Erie. 

CONTROL  OF  RIVER  ICE  FORMATION 

Diverse  structures,  such  as  ice  booms,  weirs,  and 
channel  excavations,  are  often  used  in  attempts  to 
manage  the  problems  associated  with  river  ice  forma¬ 
tion.  They  have  resulted  in  modest  success,  proving 
more  useful  for  mitigating  local  ice  problems  and  prob¬ 
lems  occurring  in  small  rivers  and  streams.  For  many 
rivers,  however,  ice  problems  occur  simultaneously  at 
several  locations  or  may  arise  at  any  location  and  may 
be  of  such  a  magnitude  as  to  make  structural  remedies 
alone  inappropriate.  Consequently,  a  question  arises  as 
to  what  approach  is  appropriate  to  mitigate  the.se  ice 
problems.  If  the  river  in  question  is  regulated,  one  ap¬ 
proach  might  be  to  control  water  discharge  and  stage 
such  that  in  conjunction  with  appropriate  structures, 
flow  conditions  are  suitable  for  ice-cover  formation  and 
maintenance. 

The  following  discussion  begins  with  a  brief  review 
of  the  various  structures  and  techniques  that  have  been 
implemented  to  promote  ice-cover  formation.  It  then 
focuses  on  the  feasibility  of  flow  regulation  to  control 
ice-cover  formation  on  the  Ohio  River  and  on  rivers 
generally. 

Control  methods 

Controlled  ice-cover  formation  usually  involves  a 
retention  structure  to  initiate  r  id  hold  an  ice  cover  and, 
on  occasion,  flow  modification  to  ensure  that  it  forms 
and  stays  in  place.  In  some  instances,  use  of  a  retention 
structure  alone  has  been  sufficient  to  develop  an  ice 
cover.  For  others,  flows  have  had  to  be  modified  such 
that  hydraulic  conditions  were  suitable  for  ice  covers  to 
form  at  retention  structures.  Flow  modification  may  en¬ 
tail  channel  alteration  and/or  regulation  of  flow  dis¬ 
charge  and  stage. 

Ice  retention  structures  serve  as  surrogate  ice  edges, 
or  barriers,  that  both  initiate  and  provide  downstream 
support  to  accumulation  ice  covers.  They  can  be  either 
temporarily  placed  across  a  channel,  such  as  ice  booms, 
or  they  may  be  constructed  as  fairly  massive,  permanent 
fixtures,  such  as  flow-regulation  dams.  The  diverse 
forms  of  retention  structures  that  have  been  used  to  date 
are  descri  bed  in  detail  by  Perham  (1983).  Michel  (1971) 
and  Ashton  (1986)  also  provide  useful  descriptions  of 
retention  structures.  The  St.  Lawrence  River,  more  than 
any  other  river,  has  been  invested  with  retention  struc¬ 
tures  for  developing  ice  covers  and  mitigating  jam 
formation. 

Floating  ice  booms  are  the  best  known  form  of  reten- 
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tion  structure.  Perham  ( 1 983)  discusses  them  at  length 
and  illustrates  their  application  at  various  rivers.  Booms 
have  several  advantages  over  other  structures: 

•  They  are  comparatively  inexpensive  to  build  and 
install, 

•  They  can  be  removed  for  ice-free  flow  conditions, 
and 

•  They  can  be  made  navigable. 

However,  they  have  to  be  used  with  due  attention  to 
flow  conditions  and,  because  of  their  light  construction, 
to  the  forces  likely  to  be  exerted  against  them. 

At  locations  w  here  severe  jams  perennially  develop 
and  where  the  limits  of  water-level  fluctuation  are 
known  and  suitably  narrow,  it  is  sometimes  advanta¬ 
geous  to  construct  fixed  booms.  /  "  advantage  of  these 
is  that  they  do  not  need  to  be  installei'.  each  winter,  but 
remain  passively  in  place.  One  such  boom  is  being  used 
to  form  and  hold  ice  covers  on  the  St.  Lawrence  River 
at  Montreal.  It  consists  of  fixed  piers  (somewhat  akin  to 
bridge  piers)  between  which  float  buoyant  beams.  Pa  riser 
et  al.  (1966)  an  1  Lawrie  (1972)  describe  its  operation. 
A  similar  installation  on  the  St.  Anne  River,  Quebec,  is 
briefly  described  by  Deck  (1984).  .A  completely  fixed 
boom  is  used  to  assist  in  ice -cover  formation  over  the 
Sigalda  Reservoir  in  Iceland  (Perham  1 983).  This  struc¬ 
ture,  however,  comprises  an  8-m-deep  continuous  beam 
that  spans  the  mouth  of  the  Tungnaa  River,  an  outlet  of 
Sigalda  Reservoir. 

Insteadof  providing  acomplete  barrier  to  form  an  ice 
cover,  it  is  sometimes  feasible  to  use  isolated  structures, 
such  as  small  artificial  islands  and  piers  or  pilings,  to 
form  ice  covers.  At  some  locations,  notably  those  where 
large  ice  floes  or  large  masses  of  frazil  slush  are  to  be 
held,  it  may  be  possible  to  precipitate  jamming  between 
structures  so  that  the  ice  acts  as  its  own  banier.  When 
frozen  into  an  ice  cover,  these  structures  also  serve  as 
point  restraints  for  the  ice  cover,  helping  to  maintain  it 
throughout  we  winter. 

The  mechanics  of  frazil  ic  ’  growth  and  transport 
have,  on  occasion,  suggested  methods  for  controlling 
ice-cover  fo.mation.  Two  that  are  still  under  investiga¬ 
tion  entail  trapping  frazil  ice  while  it  is  still  in  the  active 
phase  (water  is  still  supercooled)  and  initiating  ice- 
coverformation  close  to  the  zoneof  initial  frazil  growth. 
One  technique  (Perham  1981, 1983, 1986)  involves  the 
use  of  collector  lines  draped  in  a  flow  to  gather  frazil  ice; 
when  loaded  with  fruzil,  the  lines  become  buoyant  and 
float  to  the  flow  surface,  freeze  together,  and  thereby 
initiate  an  ice  sheet.  Another  method  (Perham  1983, 
Foltyn  1986)  entaiis  the  use  of  fence  booms,  or  frazil 
fences,  to  collect  frazil  crystals;  when  the  fence  clogs 
with  frazil,  it  retards  the  flow  and  becomes  in  effect  a 
type  of  weir.  Both  techniques  have  met  with  mixed 
results.  TTie  main  drawbacks  are  that  they  are  limited  to 


small  rivers  and  streams  and,  in  thecaseof  frazii  fences, 
may  cause  stream-bed  erosion. 

Dams  and  weirs  are  installed  as  more  or  less  perma¬ 
nent  structures  at  especially  problematical  locations. 
Their  essential  functions  are  to  reduce  flow  velocities, 
so  that  ice  covers  can  develop  over  former  high-velocity 
reaches,  and,  if  gated,  to  act  as  a  downstream  Oanier  that 
retains  ice  covers.  They  are  often  u.sed  in  conjunction 
with  floating  ice  booms  and  can  be  designed  to  different 
levels  of  sophistication.  The  simplest  and  least  expen¬ 
sive  form  are  submerged  weirs  constructed  of  loose 
boulders  and  stones;  their  suitability  is  Pmited  to  small 
rivers  and  streams.  Gated  dams  or  barrages  for  ice  con¬ 
trol  can  be  iclatively  complex  and  commensurately 
expensive  structures.  They  are  infrequently  used,  usu¬ 
ally  as  a  last  resort  at  sites  where  jam  formation  creates 
particularly  revere  problems  that  are  of  major  economic 
consequence.  One  example  of  a  gated  ice-control  dam 
is  Iroquois  Dam  located  near  the  headwaters  of  the  St. 
Lawrence  River.  Its  operation  is  described  by  Wigle  et 
al.  ( 198 1 )  and  is  discussed  later  in  this  seciton  Dams  or 
barrages  are  intended  primarily  for  flow  regulation  but 
may  also  serve  as  retention  structures;  indeed,  uninten¬ 
tionally,  they  usually  do. 

Several  concerns  may  detract  from  the  use  of  struc¬ 
tures  as  the  sole  means  of  modifying  river  flow.  To 
begin  with,  flow  conditions  have  to  be  suitable  for  ice- 
cover  formation;  usually,  this  means  that  F  should  be 
less  than  about  0. 1 5.  Suitable  flow  conditions  may  not 
be  achieved  by  structural  means.  Further,  ice  problems 
may  arise  si  multaneously  at  various  locations  or  at  any 
location  along  some  rivers.  It  may  be  both  physically 
difficult  and  uneconomical  to  place  ice  control  struc¬ 
tures  at  all  potential  problem  locations.  The  sheer  scale 
of  large  rivers  may  preclude  the  use  of  control  struc¬ 
tures:  channels  may  be  too  wide  or  too  long,  large  vol¬ 
umes  of  ice  are  produced,  associated  ice  loads  are  large, 
and  so  on.  In  addition,  installation  of  a  permanent  ice 
control  structure  may  conflict  with  use  of  a  river.  For 
example,  booms  may  pose  navigation  hazards  in  heavi¬ 
ly  trafficked  river',  and  dams  and  channel  excavations 
may  unacceptably  affect  fish  habitat. 

An  option  afforded  by  regulated  rivers  for  overcom¬ 
ing  some  of  the  foregoing  concerns  is  flow  regulation. 
There  are  a  few  ca.ses  where  this  option  has  been  utilized 
for  controlled  ice-cover  formation.  Most  have  involved 
non-navigable  rivers  and  ice-cover  formation  in  the 
vicinity  of  rapids  where  urge  quantities  of  frazil  grow. 
An  important  consideration  is  selecting  the  appropriate 
time  to  regulate  flow.  For  some  rivers,  there  may  be  few 
opportunities  to  control  ice-cover  formation.  In  this 
regard,  a  key  parameter  to  selecting  an  appropi  iate  time 
for  ice-cover  formation  is  water  temperature  in  up¬ 
stream  river  reaches  and  reservoirs. 
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Figure  4.  Controlled  ice-cover  formation  or,  the  Lule  River  in  Sweden  (Billfalk  1984 ). 


Wigle  et  al.  ( 1981 )  describe  the  measures  taken  to 
regulate  flow  from  the  outlet  of  Lake  Ontario  to  elimi¬ 
nate  excessive  frazil  ice  growth  and  to  prevent  a  frazil 
jam  from  occurring  in  the  International  Rapids  section 
of  the  St.  Lawrence  River,  Control  measures  were  im¬ 
plemented  in  order  to  prevent  flooding  and  to  mitigate 
the  deleterious  effects  of  ice  on  the  production  of  hydro- 
power  from  generation  facilities  located  immediately 
downstream  of  the  .section.  Figure  3  depicts  the  channel 
configuration.  An  ice  cover  forms  over  Lake  St.  Law¬ 


rence  as  far  as  Morrisburg.  The  control  dam  is  then  oper¬ 
ated  to  reduce  flow  (and  flow  velocities)  and  to  hold 
upstream  ice,  so  that  the  cover  can  progress  to  within  a 
few  kilometers  of  the  control  dam.  The  dam  itself,  to¬ 
gether  with  .several  ice  booms  located  upstream,  is  used 
to  establish  an  ice  cover  for  a  distance  upstream.  Shen 
and  Ho  (1986)  describe  the  use  of  a  two-dimensional 
numerical  mtxiel  to  simulate  frazil-generated  ice-cover 
formation  upstream  of  Iroquois  Dam. 

Controlled  ice-cover  formation  over  a  reach  of  the 
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Lule  River  in  northern  Sweden  is  described  by  Jensen 
(1981)  and  Billfalk  (1984).  The  river  reach  of 
is  located  on  the  lower  portion  of  the  Lule  River,  up¬ 
stream  of  one  hydropower  plant  but  downstream  of  a 
series  of  others,  as  shown  in  Figure  4.  In  this  ca.se,  flow 
regulation  for  ice  control  was  complicated  by  a  con¬ 
straint  to  maintain  power  generation.  However,  the  ap¬ 
parently  successful  approach  adopted  enables  mainte¬ 
nance  of  power  production  by  prefilling  the  lower 
reservoir  and  using  the  lower  plant  to  make  up  the  short¬ 
fall  in  power  production  from  the  upper  plants.  A 
several-day  reduction  in  flow  discharge  from  the  upper 


portion  of  the  river  and  use  of  a  boom  were  effective  in 
forming  an  ice  cover.  Earlier  attempts  at  flow  modifica¬ 
tion  by  means  of  channel  enlargement  alone  proved  to 
be  ineffective  in  enabling  ic<--covei  formation. 

Controlled  ice  covers  have  been  formed  at  one  loca¬ 
tion  on  the  Allegheny  River,  a  tributary'  to  s  -  Ohio 
River  (Fig.  5),  by  means  of  flow  regulation  and  an  ice 
boom.  Deck  and  Gooch  (1984 )  describe  the  successful 
implementation  of  an  ice  boom  to  form  and  hold  a  frazil¬ 
generated  ice  cover  over  several  kilometers  of  the  Alle¬ 
gheny  River  at  Oil  City.  The  purpose  for  forming  the  ice 
cover  was  to  suppress  e.Kcessive  frazil  ice  production 
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and  consequent  jam  occurrence  immediately  down¬ 
stream  of  Oil  City.  An  important  procedure  in  the  cover 
formation  is  regulation  of  flow  discharge  from  Kinzua 
Reservoir  upstream  of  the  boom. 

Controlled  ice-cover  formation 
for  the  upper  Ohio  River 

The  principal  impetus  for  the  present  study  is  the 
need  to  prevent  the  occurrence  of  large  ice  jams  in  stage- 
regulated  pools  of  the  upper  Ohio  River,  here  taken  to 
be  the  Ohio  River  upstream  of  Louisville,  Kentucky. 
Figure  5  illustrates  a  portion  of  the  Ohio  River  in  the 
general  vicinity  of  Pittsburgh  and  its  tributary  rivers,  the 
Allegheny  and  the  Monongahela.  In  the  past,  enormous 
damage  has  been  caused  when  large  jams  have  formed, 
led  to  flooding,  subsequently  broken  loose,  and  collided 
with  riverside  structures.  In  most  instances,  the  jams 
have  resulted  from  excessive  frazil  ice  production  and 
from  uncontrolled  ice-cover  breakup. 

A  survey  of  ice  conditions  along  the  upper  Ohio 
River,  which  was  conducted  as  part  of  the  present  study, 
leads  to  the  conclusion  that  the  most  promising  method 
of  preventing  ice  jams  involves  the  following  approach: 

Control  ice-cover  formation  through  regulation  of 
river  flow  so  that  optimal  flow  conditions  prevail  for 
rapid  formation  and  subsequent  maintenance  of  an 
accumulation  ice  cover  over  river  reaches  in  which  po¬ 
tentially  large  amounts  of  frazil  ice  may  grow.  Flow 
discharge  and  stages  would  be  regulated  by  controlling 
flow  releases  from  reservoirs  and  flow  stages  at  river 
dams.  Accumulation  covers  would  be  formed  of  frazil 
ice  pans  and  floes  and,  if  appropriate,  broken  ice  con¬ 
veyed  from  upstream.  Existing  dams,  augmented  where 
needed  by  navigable  ice  booms,  would  serve  as  reten¬ 
tion  structures  for  the  development  of  accumulation 
covers. 

This  approach  necessarily  implies  coordination  of 
flow  regulation  for  major  portions  of  the  upper  Ohio 
River  and  its  tributaries.  Before  assessing  its  feasibility 
for  pools  comprising  the  Ohio  River,  it  is  in  order  to 
present  an  overvie  w  of  flow  and  ice  conditions  along  the 
upper  Ohio. 

The  Ohio  River  and  its  two  principal  upstream  trib¬ 
utaries,  the  Allegheny  and  the  Monongahela,  are  major 
arterial  waterways  that  serve  centers  of  population  and 
industry  in  the  northeastern  United  States.  The  three 
rivers  are  heavily  used  for  shipping  and  diverse  water- 
consumption  purposes.  They  are  regulated,  at  least  par¬ 
tially,  to  meet  these  demands  as  well  as  to  prevent  flood¬ 
ing  of  adjoining  land.  Schematically,  the  rivers  can  be 
portrayed  as  comprisi  ng  series  of  navigable  pools,  punc¬ 
tuated  by  COE-controlled  dams  and  navigation  locks 
that  are  linked  to  an  array  of  flood-control  reservoirs  in 
headwater  regions  or  tributary  basins.  This  view  isillu.s- 


trated  in  Figure  6.  Tlie  pools  forming  the  Ohio  and  the 
Monongahela  rivers  are.  with  one  or  two  exceptions, 
stage-regulated  by  means  of  gated  dams.  Upstream  of 
Louisville  on  the  Ohio  River  there  are  1 4  of  these  i  nstal- 
lations.  The  Monongahela  has  nine  lock  and  dam  instal¬ 
lations.  The  Allegheny  River  is  not  stage-regulated,  but 
is  broken  into  a  series  of  eight  ungated  dams  and  locks; 
its  upper  half  is  without  dams.  I n  each  ri  ver.  the  dams  are 
intended  toprovide  navigation  channels  that  are  at  least 
about  3  m  (9  ft)  deep.  The  pools  have  negligible  storage 
capacity  and  are  not  used  for  flood  control. 

Approximately  3G%  of  flow  passing  through  the 
upper  Ohio  River  is  discharge-regulated  by  means  of  25 
flood-control  reservoirs  that  are  linked  to  one  of  the 
three  rivers;  15  reservoirs  are  operated  by  the  Corps  of 
Engineers.  The  remaining  discharge  enters  the  rivers  as 
unregulated  flows  from  numerous  small  rivers  and 
streams. 

A  somewhat  unusual  hydrologic  feature  of  the  three 
rivers  is  that  their  maximum  monthly  discharges  usual¬ 
ly  occur  during  winter  and  early  spring,  overlapping 
with  the  period  when  they  attain  minimum  water  tem¬ 
peratures  and  form  ice.  In  la.-'ge  part,  this  feature  is 
attributable  to  the  comparatively  dynamic  character  of 
winter  climate  at  the  watersheds  through  which  the 
rivers  flow.  It  is  not  uncommon  for  cold,  northerly  air 
masses  to  alternate  with  much  warmer,  moist  air  masses 
from  the  south  and  east.  The  interaction  of  these  air 
masses  is  highly  unsteady,  so  that  during  winter  the 
rivers  experience  cycles  in  prevailing  air  temperature 
and  precipitation.  Daly  (1984b)  presents  hydrologic 
and  ice  information  on  the  Ohio  River. 

Because  of  the  variable  climatic  conditions,  ice  re¬ 
gimes  on  the  rivers  are  remarkably  dynamic;  ice  forms 
during  periods  of  frigid  air,  and  ice  covers  break  up 
during  subsequent  periods  of  warm  air  accompanied  by 
rain  or  snowmelt  runoff.  Ice  conditions  on  the  Ohio 
River  have  been  observed  and  recorded  for  over  1 00 
years  by  the  U.S.  National  Weather  Service  (ORDIC 
1978).  Its  records  show  that  large  quantities  of  frazil  ice 
may  form  in  the  pools  as  well  as  border  ice  along  the 
banks.  Although  ice  forms  virtually  every  winter  on  the 
Allegheny  and  the  Monongahela  rivers  and  the  upper¬ 
most  pools  of  the  Ohio,  the  frequency  of  ice  occurrence 
diminishes  with  downstream  distance,  as  the  river  gen¬ 
erally  follows  a  southwesterly  course.  For  example,  at 
Cincinnati,  about  500  miles  downstream  of  the  conflu¬ 
ence  of  the  Allegheny  and  Monongahela  rivers,  there  is 
ice  on  the  Ohio  an  average  of  7  out  of  1 0  winters;  how¬ 
ever,  the  ice  is  usually  in  the  form  of  drifting  frazil  pans 
and  floes,  and  the  river  is  only  frozen  over  about  1 4%  of 
winters. 

Detailed  survey  charts  of  ice  conditions  on  the  Ohio, 
the  Allegheny,  and  the  Monongahela  were  prepared  by 
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Figure  6.  Hydraulic  schematic  of  the  Ohio,  Allegheny,  and  Monongahela  rivers. 


the  COE  for  the  winter  of  1 984-85  (Gatto  et  al.  1 986). 
The  charts  reveal  that,  for  the  comparatively  mild  win¬ 
ter  of  1984-85,  the  Allegheny  and  Monongahela  rivers 
and  the  Ohio  River  upstream  of  Emsworth  were  covered 
for  periods  by  fragmented  ice  covers  consisting  of  frazil 
ice  pans  and  ice  blocks.  They  also  reflect  the  rather 
dynamic  nature  of  ice  formation  on  the  Ohio,  On  aver¬ 
age,  the  Ohio  River  remains  navigable  by  tow-barge 
vessels  about  66%  of  winters.  Navigation  on  the  Al¬ 
legheny  River  usually  ceases  during  the  months  of 
January  and  February.  The  Monongahela,  located  south 
of  the  Allegheny  River  and  flowing  northward,  typical¬ 
ly  remains  navigable  throughout  winter. 

Significant  amounts  of  ice  may  enter  the  upper 
reaches  of  the  Ohio  River  and  all  of  the  Allegheny  and 
Monongahela  rivers  from  their  numerous  small  tribu¬ 
taries.  Being  relatively  small  water  bodies,  such  tribu¬ 
taries  cool  rapidly  and  produce  ice  in  most  winters.  The 
ice  on  them  also  breaks  up  and  is  conveyed  to  the  Ohio, 
often  before  the  ice  on  that  river  has  broken  up. 


Prominent  among  the  problems  arising  from  ice 
formation  on  the  Ohio,  Allegheny,  and  Monongahela 
rivers  are  those  attendant  to  ice  jams;  a  summary  of  ice 
problems  on  these  rivers  is  given  by  Zufelt  and  Calkins 
(1985).  For  the  Ohio  River,  the  economic  consequences 
of  these  problems  are  exacerbated  by  its  heavy  use  and, 
relatedly,  by  its  passage  through  extensively  inhabited 
and  industrialized  regions.  The  Markland  incident,  re¬ 
counted  in  detail  by  the  Ohio  River  Division  Ice  Com¬ 
mittee  (ORDIC  1978),  provides  stark  testimony  to  the 
damage  that  can  be  caused  by  a  large  Jam  in  the  Ohio 
River.  The  effects  of  dynamic  climate  and  the  lack  of 
synchronized  regulation  of  major  segments  of  these 
rivers  have  led  to  the  troublesome  occurrence  of  ice 
jams.  Some  jams  recur  at  several  locations  on  the  Ohio 
and  the  Allegheny.  As  indicated  in  Figure  6,  these  loca¬ 
tions  often  coincide  with  the  confluence  of  tributary 
rivers  or,  as  noted  by  Zufelt  and  Calkins  (1985),  with 
bends,  bridge  crossings,  islands,  and  channel  constric¬ 
tions. 
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Each  lock  and  dam  installation  is  a  potential  jam  site; 
viz.  the  Markland  incident.  If  there  is  no  coordination 
between  flow  and  ice  releases  along  a  waterway  com¬ 
prising  a  series  of  lock  and  dam  installations  and  reser¬ 
voirs,  there  is  the  potential  for  a  jam  to  occur  at  any  in¬ 
stallation. 

NUMEWCAL  SIMULATION  OF 
ICE-COVER  FORMATION 

With  appropriate  regulation  of  flow,  lock  and  dam 
installations  could  be  used  as  retention  structures  from 
which  accumulation  ice  covers  could  be  formed.  We 
will  now  examine  this  concept  further  through  numer¬ 
ical  simulation  of  accumulation  cover  formation. 

Numerical  simulation  of  accumulation  ice-cover  for¬ 
mation  resulting  from  frazil  ice  generation  involves 
computation  of  the  following: 

•  Flow  profile, 

•  Water  temperature, 

•  Rate  of  frazil  ice  growth, 

•  Rate  of  ice  cover  progression  upstream,  and 

•  Ice  cover  thickening  by  accumulation  and  static 
(thermal)  ice  growth. 

In  thepresentstudy,  one-dimensional  unsteady-flow 
equations  are  used  to  determine  the  first  three  items.  Ice- 
cover  thickness  profiles  are  calculated  using  algebraic 
relations  involving  a  stability  criterion  for  ice  pans,  or 
floes,  arriving  at  the  leading  edge  of  an  ice  cover. 

An  outline  is  presented  here  of  the  formulations, 
boundary  and  initial  conditions,  and  computational 
schemes  used  in  developing  the  numerical  simulation. 
The  following  general  assumptions  are  made  about 
frazil  ice  growth  and  accumulation  ice-cover  forma¬ 
tion: 

(i)  Flow,  frazil  ice  growth,  and  ice-cover  progres¬ 
sion  can  be  described  in  terms  of  depth-averaged, 
one-dimensional  formulations; 

(ii)  The  two-layer  hypothesis*  is  appropriate  for  de¬ 
scribing  ice-covered  river  flow; 

(iii)  River  discharge  is  kept  constant  during  ice-cover 
formation; 

(iv)  Ice  covers  form  at,  and  develop  upstream  from, 
stage-regulation  dams  where  water  stage  is  held 
constant; 

(v)  Frazil  ice  forms  in  the  reach,  or  pool,  of  a  river 
over  which  an  ice  cover  is  to  be  formed;  no  frazil 
or  broken  ice  is  transported  from  the  reach  above 
an  upstream  dam  (which  in  effect  could  be  con¬ 
sidered  as  also  becoming  ice  covered);  and. 


•  The  iwo-layer  hypothesis  holds  that  ice-covered  flows  in  formerly 
open  channels  can  be  treated  as  two  free-surface  flows:  one  with  the 
channel  bed  as  its  base,  the  other  with  the  underside  of  the  ice  cover 
as  its  base. 


(vi)  Border  ice  growth  does  not  significantly  narrow 
flow  surface  during  accumulation  ice-cover  for¬ 
mation. 

Other  assumptions  are  noted  within  the  description 
of  each  component  of  the  simulation  model. 

One  issue  that  is  a  major  stumbling  block  in  comput¬ 
ing  ice-cover  formation  from  frazil  ice  is  the  lack  of  in¬ 
formation  concerning  lengths  of  flow,  or  flow  periods, 
required  for  frazil  in  the  form  of  individual  crystals  to 
agglomerate  as  surface-floating  frazil  pans.  Some  re¬ 
search  effon  has  been  invested  in  examining  the  “rise 
velocity”  of  frazil  crystals  and  floes  (Park  and  Gerard 
1984,  Wuebben  1984),  but  the  results  do  not  fully 
address  this  issue.  Another  issue  concerns  the  propor¬ 
tions  of  frazil  ice  grown  that  are  transformed  to  frazil 
pans  or  frazil  slush.  Here,  toi,  there  are  few  quantitative 
guidelines  lobe  implemented  in  a  numerical  simulation 
model.  To  circumvent  these  issues,  while  remaining 
cognizant  of  their  significance,  two  further  assumptions 
were  made  concerning  accumulation  cover  formation. 
First,  it  is  assumed  that  all  frazil  ice  grown  upstream  of 
an  advancing  cover  goes  to  form  that  cover.  Second,  as 
a  means  of  accounting  for  mass  of  frazil  ice  grown, 
cover  progression  was  calculated  such  that  a  cover 
advances  in  incremental  steps  equivalent  to  the  compu¬ 
tational  step  length.  Ax.  For  each  incremental  length  of 
advance,  the  cover  thickens  until  it  attains  a  thickness 
equivalent  to  that  of  a  typical  pan,  which  is  assumed  to 
be  102  mm  (4  in.).  Once  the  incremental  advance  has 
reached  this  thickness,  the  stability  criterion  associated 
with  the  juxtaposition  of  frazil  ice  pans  is  applied.  An¬ 
other  way  of  viewing  the  simulation  is  that  frazil  pans 
were  assumed  to  accumulate  at  the  head  of  an  advancing 
cover  in  such  a  way  that  the  head  advances  a  computa¬ 
tional  step.  Ax,  before  the  stability  criterion  is  applied 
for  accumulation  cover  formation  by  juxtaposition.  If 
this  criterion  is  not  met,  the  cover  thickens,  and  flow 
conditions  are  modified  until  it  is  met.  Admittedly,  the 
simulation  is  not  precise,  but  within  the  context  of  a  pre¬ 
liminary  feasibility  study  it  is  tenable  and  yields  mean¬ 
ingful  information. 

Flow  profile 

Governing  equations 

Profiles  of  flow  through  rivers  with  floating  ice 
covers  (Fig.  7)  are  described  here  using  one-dimension¬ 
al  equations  of  gradually  varied  unsteady  flow.  These 
equations  are: 

Conservation  of  liquid-water  mass, 

^-1-— =  0,and  (1) 

dx  dt 
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Figure  7.  Definition  sketch. 


Conservation  of  flow  momentum. 


+ 

dx\Al 


gA^+gAS(=Q 

ax 


Initial  and  boundary  conditions 

The  upstream  boundary  condition  is  constant  water 
discharge  Q: 

(2(0.0  =  (2„  (5) 


where  Q  =  water  discharge 

A  =  cross-sectional  area  of  flow 
S{  =  friction  slope 
g  =  gravitational  acceleration 
X  -  streamwise  distance 
t  =  time  _ 

H  =  z  +  d  +  r\ 
z  -  bed  elevation 
d  =  flow  depth 

T]  =  jjT]  =  equivalent  thickness  of  ice  cover 
Sj  =  specific  gravity  of  ice 
q  =  ice-cover  thickness. 

Friction  slope,  Sf,  can  be  expressed  in  terms  of 
Manning’s  coefficient  as 

V  =  Jlr^S(^  (3) 

nc 

where  V  =  mean  velocity  of  flow, 

R  -  hydraulic  radius,  and 

«£  =  composite  Manning’s  coefficient. 

In  this  report,  for  ice-covered  flow,  the  composite  Man¬ 
ning’s  coefficient  is  taken  to  be 

nM0Mn>'^  +  n^'^)r  (4) 

in  which  nj,  and  n\  are  the  Manning’s  coefficients 
related  to  flow  layers  associated  with  the  channel  bed 
and  ice  cover,  respectively  (Uzuner  1975).  Note  that 
loss  of  liquid  water  by  freezing  to  ice  is  neglected  in  eq 
1 ,  as  the  volumetric  rate  of  ice  growth  is  much  smaller 
than  the  water  discharge. 


The  downstream  boundary  condition  is  constant 
water  depth  at  the  downstream  control  section, 

d{e,t)  +  qff  ,0  =  Yt  (6) 

in  which  t  is  the  length  of  river  reach.  The  initial  flow 
profile  is  obtained  by  a  stabilization  process  in  which 
the  unsteady  flow  equations  for  fixed  boundary  condi¬ 
tions  are  solved  until  an  assumed  flow  profile  becomes 
steady. 

Numerical  scheme 

A  weighted,  implicit  finite-difference  scheme  (Cunge 
et  al.  1980),  first  developed  by  Preissmann  in  1960,  is 
used  to  solve  simultaneously  eq  1  and  2.  In  accordance 
with  Preissmann’s  computational  scheme  (Fig.  8),  on 


Figures.  A  weighted,  implicit  finite-difference  scheme. 
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omitting  the  superscript  J  for  convenience  with  the 
understanding  that^i’jt.ri  without  a  superscript  corre¬ 
sponds  to  the  value  of  f(x,t)  at  the yth  time  step,  any 
variabieyi'x.ri  and  its  derivatives  are  discretized  as 

/U)  =  p/i^,  +  ( 1  _  p)/i  +  e  [pA/.+i  -H  { 1  -  p)  A  /i]  (7) 

s  J-  1/i.i  -  /i  +  e (A/,>,  -  A/.)]  (8) 

ox  Ax 

^  =  X[pA/i.,  +  (l-p)A^]  (9) 

dt  At 

where  A  /j  =  -fJ 

Ax  =  length  increment 
At  =  time  increment 
(  =  node  along  x-axis 
j  =  node  along  r-axis 
6  =  temporal  weighting  factor 
P  =  spatial  weighting  factor. 

Cunge  and  Perdreau  (1973)  suggested  a  practical  range 
for  0  of  0.6  to  1 .0.  In  the  present  study,  we  use  0  =  0.65 
and  P  =  0.5;  using  slightly  different  values  does  not 
affect  the  results  of  the  present  study. 

Equations  1  and  2,  on  substituting  eq  7  through  9,  can 
be  written  symbolically 

OiAyi+i  +  biAQi^i  =  aAyi  +  diAQi  +  ei  (10) 

OfAyi+i  +  b,'Aj2t+i  =  C(Ayi  +  d{AQi  +  (11) 

in  which  y  =  z  +  d. 

The  coefficients  aj,  bj,  cj,  dj,  ej,  Of,  b(,  C(,  d(,  and  «,-are 
given  in  Appendix  A.  When  eq  8  and  9  are  applied  to  the 
first  (.N~  1 )  grid  points  along  the  length  of  the  channel  {N 
is  the  total  number  of  grid  points),  a  system  of  two  (N- 
1)  linear  equations  involving  2N  unknowns,  AQ;  and 
Ayj,  results.  With  two  boundary  conditions  (one  at  the 
upstream  boundary  and  the  other  at  the  downstream 
boundary  for  subcritical  flows),  the  number  of  equa¬ 
tions  matches  the  number  of  unknowns,  and  hence  the 
system  of  algebraic  equations  can  be  solved  to  obtain 
the  solution  (i.e.,  y/^*,  fif"* )  at  the  next  time  step.  The 
resulting  system  of  equations  is  solved  using  the  dou¬ 
ble-sweep  method.  It  is  assumed  that  q  is  constant 
within  At. 

Water-temperature  variation 

Governing  equation 

Streamwise  variation  of  depth-averaged  water  tem¬ 
perature  can  be  expressed  as  a  one-dimensional  advec- 
tion-diffusion  equation  (Brocardand  Harlemann  1976) 
in  terms  of  conservation  of  thermal  energy; 


—  (pCp4  r)H - {QpCpT) 

dt  dx 


a  /  ^7* 

=  -^/V£xpCp—  +fl0wa  (12) 

3x\  dx, 

where  T  =  depth-averaged  temperature  of  water 
B  =  channel  width 
p  =  density  of  water 
Cp  =  specific  heat  of  water 

=  longitudinal  dispersion  coefficient 
=  net  heat  flux  from  water  to  air,  per  unit  sur¬ 
face  area  of  flow. 

For  Q= constant,  and  neglecting  longitudinal  disper¬ 
sion,  eq  1 2  can  be  reduced  to  a  simplified  form; 

—  9wa  (]3) 

3r  3x  pcprf 

Net  heat  flux  <|>^a  is  taken  to  be  a  linear  function  of  the 
air  and  water  temperature  difference,  such  that 

(|)wa  — /*wa  (Ta  —  7^  (14) 

where  is  the  heat  transfer  coefficient  at  the  water/air 
interface  and  is  the  air  temperature.  Herein,  = 20 
'W/m'^J°C  is  used,  following  the  recommendations  of 
Ashton  (1986)  and  Shen  and  Ho  ( 1 986).  Substitution  for 
(l),^,a  from  eq  14  into  eq  13  yields 

^+V—  =  Jb!t^in-T)  (15) 

dt  dx  pcvd 

Initial  and  boundary  conditions 

The  following  quantities  are  taken  to  be  constant; 
initial  water  temperature  along  the  river  reach  of  inter¬ 
est,  temperature  of  in-flow  water  with  time,  and  air 
temperature  both  above  the  river  reach  and  with  time.  In 
other  words. 


nx,o)=To 

(16) 

no,t)  =  T„ 

(17) 

'P 

It 

o 

(18) 

Numerical  scheme 

The  method  of  characteristics,  with  Holly-Preiss- 
mann's  4th-order  interpolation  scheme  (Holly  and  Pre- 
issmann  1977),  is  used  to  solve  eq  15,  which  can  be 

transformed  into  a 

set  of  two  ordinary  differential 

equations; 
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X 

Figure  9.  A  computational  scheme formethod  of  charac¬ 
teristics. 

QL+a(T-n]  =  0  (19) 

Dt  ^  ’ 

along  a  characteristic  curve  (Fig.  9)  defined  as 

dx.  =  v  (20) 

dt 

In  eq  1 9,  a  =  h^J(pCpd).  Integration  of  eq  1 9,  from 
ti  to  along  a  characteristic  curve  gives  (Fig.9) 

a{T-T^)dt  =  0.  (21) 

Using  trapezoidal  integration,  eq  21  can  be  approxi¬ 

mated  as 

7/^'  =  ^[(ar)f'  -(arjf  +  {o.T)^-{oim] 

(22) 

Unknown  values  are  evaluated  using  eq  22,  and 
values  of  7^  are  interpolated  using  Holly-Preissmann’s 
4th-order  interpolation  scheme,  which  introduces  neg¬ 
ligible  numerical  diffusion  and  dispersion. 

Frazil  ice  growth 

Governing  equation 

It  is  assumed  that  the  thermal-energy  loss  to  air  of 
freezing  water  is  balanced  by  the  thermal  energy  gener¬ 
ated  by  frazil  ice  growth.  The  thermal  energy  generated 
per  unit  volume  of  frazil-laden  water  is  PiZ-jC,  with  pj  = 
density  of  ice,  Lj  =  latent  heat  of  water  fusion,  and  C  = 
volumetric  concentration  of  the  frazil.  By  replacing 
pCp7  with  -PiZ-jC,  eq  15  can  be  used  to  determine  con¬ 
centration  of  frazil  ice  grown; 

(23) 

dt  dx  fuLid^  ' 
in  which  7f  is  the  freezing  temperature  of  water,  herein 
taken  as  0°C. 


Initial  and  boundary  conditions 

As  mentioned  above,  it  is  assumed  that  there  is  zero 
influx  of  frazil  ice  entering  the  upstream  boundary  of 
the  river  reach  under  consideration: 

C(0,f)=:0  (24) 

Because  initial  water  temperature  in  an  entire  river 
reach  is  assumed  to  be  above  0°C,  initial  concentration 
of  frazil  ice  is  taken  as  being  zero; 

C[x,0)=0  (25) 

Numerical  scheme 

Equation  23  is  solved  for  C  using  the  same  method 
of  characteristics  used  to  solve  for  7.  Integration  of  eq 
23  along  the  characteristic  curve  (eq  20)  from  0  to 
(Fig.  9)  yields 


[{erf)>'-(e7’JA'  +  (e7f)5-(er^4]  (26) 

in  which  e  =  h^Jp^l^d.  Values  of  are  interpolated 
using  Holly-Preissmann’s  4th-order  interpolation 
scheme. 

Ice-cover  progression 

Itisassumed  that,  once  initiated,  an  icecoverprogress- 
es  upstream  in  incremental  steps  of  length  Ajc  through 
accumulation  of  incoming  frazil  ice.  When  Froude 
number  F  at  the  leading  edge  of  the  ice  cover  is  less  than 
a  critical  Froude  number  F^,  for  the  formation  of  ice 
cover  by  juxtaposition  of  frazil  ice  pans,  the  leading 
edge  progresses  upstream  as  a  single  layer  of  pans. 
Based  on  field  observations  reported  by  Shen  et  al. 
(1984),  F^  is  taken  here  to  be  0.05.  The  minimum  thick¬ 
ness  of  accumulation  cover  so  formed  is  the  thickness  of 
individual  ice  pans,  /ij,  here  assumed  to  be  102  mm  (4 
in.). 

When  Fexceeds  F^,  water  flow  submerges  incoming 
frazil,  which  is  transported  beneath  and  deposits  along 
the  underside  of  the  ice  cover.  This  leads  to  a  gradual  in¬ 
crease  of  ice  cover  thickness  in  the  vicinity  of  the  ice 
front.  The  velocity  that  is  required  for  ice  pans  or  floes 
to  submerge  below  the  ice  cover  increases  with  the  in¬ 
creasing  thickness  of  the  ice  cover.  For  a  given  flow 
depth  and  velocity  upstream  of  an  ice  front,  a  limiting 
value  exists  for  the  accumulation  thickness  of  an  ice 
cover  during  its  progression.  Provided  F  does  not  ex¬ 
ceed  a  maximum  value  F^^^,  such  that  pans  do  not  de¬ 
posit  and  thicken  an  ice  cover,  an  ice  covercan  progress 
upstream.  The  resulting  local  thickness,  rj,  can  be  com¬ 
puted  from  an  expression  proposed  by  Michel  ( 1 978)  as 
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Shen  and  Chiang  1984).  By  equating  heat  conduction 
(27)  through  an  ice  sheet  to  heat  efflux  from  its  top  surface, 
7s  can  be  determined  as 


where  =  tp  +  (1  -  ep)e  =  overall  porosity  of  the  ice 
accumulation, 

Cp  =  porosity  of  the  ice  accumulation  between 
ice  pans,  and 

e  =  porosity  of  individual  ice  floes. 


A  value  of  0.09  is  used  (following  Shen  et  al.  1984)  for 
the  average  value  of  which  corresponds  to  the 
maximum  value  of  x\ld. 

The  change  in  ice  thickness,  AtIj,  of  a  reach  Ax  long 
in  period  Ar  is 


Ati,= 


gCAr 

BAx 


(28) 


If  F  exceeds  a  maximum  value  frazil  passes  be¬ 

neath  the  ice  cover  and  is  swept  downstream  until  flow 
velocity  becomes  low  enough  for  the  frazil  to  rise  and 
come  to  rest  beneath  the  ice  cover,  possibly  forming  a 
hanging  dam  of  frazil. 


Thermal  growth  of  ice  cover 

Governing  equation 

As  an  accumulation  cover  progresses  upstream,  it 
solidifies  as  porous  frazil  ice  is  thermally  transformed  to 
monolithic  ice  and  thermally  thickens.  Thermal  thick¬ 
ening  is  simulated  using  the  following  thermal  energy 
equation: 


Pi  Li  ^  ^  {7f  -  fs)  -  Awi  {7  -  7f)  (29) 

where  7^  =  temperature  on  the  top  surface  of  the  ice 
cover, 

k  =  thermal  conductivity  of  ice,  and 
h^^  =  heat  transfer  coefficient  at  the  water/ice 
interface. 

The  first  and  second  terms  on  the  righthand  side  of  eq 
29  are  ice-sheet  conduction  of  thermal  energy  and 
thermal  energy  transfer  to  the  underside  of  the  ice-sheet, 
respectively. 


j’  _ ^7f'4*  fljaTjl  (31) 

*  hafX  +  k 

NUMERICAL  RESULTS 

The  simulation  model  was  exercised  to  determine 
the  times  and  corresponding  volumes  of  ice  required  to 
form  frazil-generated  accumulation  covers  over  regu¬ 
lated  river  pools.  The  parameters  varied  were  pool 
length,  pool  bottom  slope,  and  water  discharge. 

The  time  required  to  form  an  ice  cover  by  frazil  ice 
accumulation  over  a  certain  length  of  a  river  in  a  naviga¬ 
tion  pool  depends  on  a  number  of  variables  (Table  1 ) 
and  coefficients  (Table  2)  as  well  as  ice  and  water  prop¬ 
erties  (Table  3).  The  numerical  experiments  were  car¬ 
ried  out  with  fixed  values  of  the  coefficients  and  the  ice 
and  water  properties,  which  are  also  included  in  Tables 
2  and  3.  Though  the  computer  program  (Appendix  B) 
was  developed  for  general  longitudinal  and  cross-sec¬ 
tional  channel  shapes,  each  stage-regulated  pool  is 

Table  1.  Geometric  and  flow  variables. 


Variable 

Hannibal 

pool 

Montgomery 

pool 

Pool  length,  t  (km) 

66 

30 

Channel  width,  B  (m) 

335 

305 

Downstream  flow  depth,  Y/  (m) 

ll.S 

8.7 

Composite  Manning’s  coefficient, 

0.033 

0.033 

Channel-bottom  slope. 

1.17  xl(H 

1.97  X  10^ 

Table  2.  Values  of  empirical  coefficients. 

Assigned 

Empirical  coefficient _ Symbol  value 


Critical  Froude  number  for  cover  formation 


by  juxtaposition  of  ice  pans 

Maximum  Froude  number  for  ice-cover 

0.05 

progression 

0.09 

Heat-transfer  coefficient  at  air/ice  interface 

20  W/m2/°C 

Heat-transfer  coefficient  at  water/air  interface 

^wa 

20  W/m2/°C 

Thermal  conductivity  of  ice 

k 

2.22  W/mTC 

Latent  heat  of  fusion  of  ice 

k 

333  J/g 

Numerical  scheme 

Equation  29  is  used  to  determine  ice  cover  thicken¬ 
ing.  Assuming  a  quasi-steady  state  during  a  time  step  Ai, 
change  in  ice  thickness  (Aq)  is 

Aq,  =  -M.  [h;^{T,-  7a)-  Awi  (7-  7f)]  (30) 

pi  Li 

where  h,g  is  the  heat  transfer  coefficient  at  the  air/ice 
interface,  here  taken  to  be  20  W/m^FC  (Ashton,  1986, 


Table  3.  Ice  and  water  properties. 


Physical  properly 

Symbol 

Assigned 

value 

Specific  heat  of  water 

Cp 

4.2  kJAgTC 

Porosity  of  individual  ice  pans 

e 

0.4 

Porosity  of  ice  accumulation 

0.4 

Specific  gravity  of  ice 

0.917 

Density  of  water 

p 

10^  kg/m^ 

Density  of  ice 

Pi 

917  kg/m3 
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assumed  to  have  constant  values  of  longitudinal  slope 
and  Manning’ s  coefficient  and  to  be  of  constant  rectan¬ 
gular  cross-section. 

niiistration  of  simulated  ice-cover  formation 

As  an  illustration  of  what  the  numerical  model  sim¬ 
ulates.  Figures  10  and  1 1  show,  respectively,  predicted 
profiles  of  ice  cover  and  flow-depth  profile  in  a  repre¬ 
sentative  regulated  pool  of  the  Ohio  River;  the  pool 
length  is  24  kra.  For  the  sake  of  clarity,  the  effect  of 
thermal  growth  of  the  ice  cover  is  not  included  in  these 
two  figures.  Initial  flow  profile,  indicated  by  the  line  (f 
= 0),  varies  from  8.5  m  (28  ft)  at  the  downstream  end  of 
the  pool  to2.8  m  (9  ft)  at  its  upstream  end.  The  unit  water 
discharge  is  2.8  m^/s/m  (30  ftVs/ft),  which  yields  F  = 
0.1 2  in  the  upstream,  uniform-flow  reach  and  F<  0.05 
near  the  downstream  end  of  the  pool,  from  about  jc  =  2 1 
to  24  km.  Frazil  ice  pans  start  accumulating  at  the 
downstream  end  of  the  pool  where  a  navigation  dam 
serves  to  retain  them  and  to  regulate  flow  stage.  The 


thickness  of  the  ice  cover  in  the  reach  from  a:  =  2 1  to  24 
km  is  102  mm  (4  in.)  during  the  early  period  of  ice 
accumulation,  because  F  within  this  reach  is  less  than 
0.05  (the  criterion  for  pan  stability  and  upstream  pro¬ 
gression  of  ice  cover)  and  because  the  ice  cover  thick¬ 
ness  is  controlled  by  simple  juxtaposition  of  frazil  pans. 

The  ice  cover  did  not  progress  continuously.  It  ad¬ 
vanced  until  time  t2,  stopped  advancing  between  times 
t2  and  /3,  resumed  upward  progression  until  lime 
when  it  stopped  again,  and  so  on,  as  shown  in  Figure  1 0. 
The  inability  of  the  ice  cover  to  progress  upstream 
between  times  /2  and  and  and  ty  is  attributable  to 
F  exceeding  0.05  at  the  leading  edge  of  the  ice  cover 
during  these  periods.  However,  incoming  frazil  was 
deposited  under  the  ice  cover,  which  caused  it  to  thicken 
and  the  resultant  water  depths  upstream  from  it  to  in¬ 
crease,  as  shown  in  Figure  1 1 .  When  Freduced  to  0.05, 
the  icecover  resumed  its  upstream  progression.  Through¬ 
out  the  formation  of  the  ice  cover,  flow  stage  is  main¬ 
tained  constant  at  the  poofs  downstream  end.  However, 


water  depth  at  the  leading  edge  of  the  ice  cover  is 
constant  at  about  4.6  m,  as  can  be  seen  in  Figure  11, 
corresponding  to  F= 0.09.  The  incoming  ice  floes  could 
no  longer  be  deposited  in  this  region  and  were  carried 
downstream  before  being  deposited  beneath  the  ice 
cover. 

It  is  evident  from  this  example  that  the  stability  cri¬ 
terion  for  cover  progression,  F  <  0.05,  significantly 
affects  the  overall  period  of  time  required  for  an  accu¬ 
mulation  cover  to  progress  upstream.  Although  the  sen¬ 
sitivity  to  this  criterion  of  rate  of  ice-cover  progression 
is  not  quantitatively  examined  in  this  report,  it  can  be 
demonstrated  that  extension  of  the  criterion  to  F<  0. 1 5, 
as  proposed  by  Ashton  ( 1 986),  would  lead  to  more  rapid 
ice-cover  progression  and  less  ice  grown. 

In  the  following  sections,  attention  is  concentrated 
primarily  on  the  overall  time  required  for  an  accumula¬ 
tion  cover  to  form  over  pools  or  reaches  of  rivers,  the 
total  volume  of  ice  grown,  and  the  thickness  profile  of 
the  resulting  accumulation  cover. 

Ice-cover  formation  in  the 
Hannibal  and  Montgomery  pools 

Numerical  experiments  were  conducted  to  deter¬ 
mine  the  time  required  for  accumulation  covers  to  form 
for  two  pools  in  the  Ohio  River,  those  immediately 
upstream  of  the  lock  and  dam  installations  at  Hannibal 
and  Montgomery  (Fig.  5).  The  idealized  dimensions 
and  geometries  of  these  pools  are  included  in  Table  1. 
The  main  difference  in  their  geometries  is  in  their  bot¬ 
tom  slope,  downstream  flow  depth,  and  the  pool  length. 
The  bottom  slope,  downstream  flow  depth,  and  length 
of  the  Montgomery  pool  are  respectively  steeper,  high¬ 
er,  and  shorter  than  those  of  the  Hannibal  pool.  Ice- 
cover  formation  was  studied  for  three  air  temperatures 


(-10,-15,  and  -20“C)  and  a  range  of  river  flows. 

Temporal  evolutions  of  accumulation  ice  covers 
over  the  Hannibal  pool,  under  an  air  temperature  of 
-20°C,  are  shown  in  Figures  12,  13,  and  14  for  three 
river  discharges:  100, 600,  and  900  m-^/s,  respectively. 
For  the  two  lower  discharges  (100  and  600  m^/s),  ice- 
cover  thickness  is  approximately  uniform  along  the 
entire  pool  length,  being  slightly  thickerthan  1 02  mm  as 
given  by  the  condition  of  pan  juxtaposition.  The  in¬ 
crease  in  thickness  above  1 02  mm  is  due  to  thermal  ice 
growth.  In  effect,  for  these  two  discharges.  F  at  the 
leading  edge  of  the  ice  cover  never  exceeded  the  as¬ 
sumed  critical  value  of  0.05. 

The  resulting  ice  cover  is  not  uniformly  thick  when 
river  discharge  is  900  mVs.  Instead  of  forming  as  a 
simple  juxtaposition  of  frazil  pans,  it  progressed  with  a 
thickened  leading  edge  because  there  falways  exceed¬ 
ed  0.05.  The  additional  thickening  of  the  ice  cover  due 
to  submergence  and  deposition  of  frazil  ice  delays  ice- 
cover  progression. 

The  time,  required  for  an  ice  cover  to  reach  the 
upstream  end  of  the  Hannibal  pool  is  a  function  of  river 
discharge  and  air  temperature,  as  shown  in  Figure  15. 
Tire  interesting  aspect  of  this  figure  is  that,  for  a  given 
air  temperature,  it  indicates  the  existence  of  an  optimum 
river  discharge  corresponding  to  a  minimum  time  re¬ 
quired  for  the  cover  to  form  over  the  pool.  Herein  such 
optimal  discharges  are  designated  as  Q^,. 

The  physical  explanation  for  the  occurrence  of  such 
optimal  discharges  is  as  follows.  The  rate  of  ice  accu¬ 
mulation  increases  with  increasing  river  discharge,  as 
can  be  seen  by  comparing  the  areas  under  a  curve  for  t 
=  constant,  say  equal  to  50  hours,  in  Figures  12  to  14. 
Thermal  growth  of  the  ice  cover  does  not  alter  the  rela¬ 
tionship  between  the  rates  of  ice  accumulation  and  river 
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Figure  12.  Ice  cover  development  over  the  Hannibal  pool  for  Q  =  WO  m^/s. 
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Figure  13.  Ice-cover  development  over  the  Hannibal  pool  for  0  =  600  ni^/s. 
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Figure  1 4,  Ice-cover  development  over  the  Hannibal  pool  for  Q  =  900  m^/s. 
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Figure  15.  Time  for  ice-cover  formation  over  the  Hannibal  pool. 
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Figure  16.  Volume  of  ice  in  ice  cover  over  the  Hannibal  pool. 


Figure  1 7.  Time  for  ice-cover  formation  over  the  Mortgomeiy  pool. 
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Figure  18.  Volume  of  ice  in  ice  cover  over  the  Montgomery  pool. 
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discharge,  though  the  respective  areas  do  show  addi¬ 
tional  thickening  that  is  attributable  to  thermal  thicken¬ 
ing.  One  would  expect  the  requisite  time  to  decrease 
with  increasing  rate  of  ice  accumulation  and,  conse¬ 
quently,  with  an  increasing  river  discharge.  The  infer¬ 
ence  is  correct  as  long  as  ice  thickness  does  not  change 
with  river  discharge,  as  would  be  the  case  when  fat  the 
leading  edge  of  the  ice  cover  never  exceeds  0.05.  This 
condition  would  be  satisfied  for  river  discharges  less 
than  a  certain  critical  value.  For  river  discharges  larger 
than  Q^,  F  at  river  reaches  toward  the  upstream  end  of 
the  pool  would  exceed  0.05,  causing  delay  in  ice  cover 
progression,  nS  explained  earlier,  and  extending  r^. 

Figure  15  also  reflects  the  fact  that  the  rate  of  frazil 
ice  growth  increases  with  decreasing  air  temperature 
and,  therefore,  decreases  with  decreasi  ng  air  temper¬ 
ature  for  a  given  river  discharge.  Interestingly,  is  not 
sensitive  to  air  temperature  because  it  is  controlled  by 
the  critical  value  of  f  =  ^  =  Q/ffyi  5^0.5)  =  o.05 

The  volumes,  V,  of  ice  contained  in  an  ice  cover 
(both  as  accumulated  ice  and  thermal  ice  growth)  over 
the  Hannibal  pool  at  time  for  various  river  discharges 
and  the  three  air  temperatures  of  -10,  -15,  and  -20“C 
are  shown  in  Figure  16.  Volumes  of  ice  grown  are 
almost  constant  tor  river  discharges  less  than  j2o  *5e- 
cause  ice  covers  form  as  a  single  layer  of  juxtaposed 
pans  with  uniform  initial  thickness.  For  discharges 
higher  than  Q^,  nonuniform  thick  accumulation  covers 
forna  such  that  V  increases  with  increasing  river  dis¬ 
charge.  The  volume  of  ice  at  time  for  a  given  dis¬ 
charge  is  not  sensitive  to  air  temperature,  because  ice- 
cover  profiles  are  primarily  controlled  by  the  flow  con¬ 
ditions.  With  decreasing  air  ternperatu'^es,  the  rate  of 
thermal  ice  growth  increases  and  time  decreases.  The 
net  effect  of  air  temperature  on  ice  volume  resulting 
from  thermal  thickening  of  the  ice  cover  is  insignificant 
during  ice-cover  progression. 

The  numerical  experiments  produce  results  for  the 
Montgomery  pool  similarto those  forthe  Hannibal  pool 
(Fig.  1 7  and  1 8).  However,  because  the  Montgomery 
pool  is  much  shorter  than  the  Hannibal  pool,  lesser 
times  are  required  to  form  ice  covers,  and  the  resulting 
volumes  of  ice  are  less. 

The  oldest  gauging  station  nearest  to  the  Montgom¬ 
ery  and  Hannibal  pools  is  located  at  Sewickly,  Pennsyl¬ 
vania  (Fig.  5).  The  average  monthly  discharge  at  this 
sta-tion  from  December  to  February  ranges  from  about 
250  to  1,250  m3/s  (8,81 1  to  44,057  ftVs).  In  compari¬ 
son,  the  numerical  simulation  indicates  that  values  of 
Qo  for  the  Hannibal  and  Montgomery  pools  are  6(X)  and 
300  mVs  (21,147  and  10,574  ft^/s),  respectively,  when 
flow  stage  is  held  at  normal  operating  level  at  each 
installation.  For  higher  stages,  large  values  of  Qq  could 
be  used  (see  Generalized  Results,  below).  Because 


for  each  pool  falls  within  the  range  of  observed  mean 
monthly  flows,  it  appears  that  flow  (primarily  dis¬ 
charge)  regulation  can  be  used  effectively  to  minimize 
the  duration  of  ice-cover  formation  ever  these  pools. 
The  duration  of  regulation  associated  with  for  each 
'■ise  is  about  70  to  100  hours  (nominally  three  to  four 
days).  An  issue  arises,  however,  for  a  river  comprising 
several  pools  of  varying  geometry,  Qg,  for  one  pool,  or 
reach,  may  not  be  the  same  as  for  others  or  the  river  as 
a  whole. 

Generalized  results 

Numerical  simulations  of  ice-cover  formation  over 
the  Hannibai  and  Montgomery  pools  indicate  that  there 
existoptimumriverdischarges,  00,  associated  with  ice- 
cover  formation  in  minimal  periods.  The  simulations 
were  extended  to  determine  values  of  Qq  for  a  range  of 
regulated  river  pools  under  the  following  conditions: 

(i)  Air  temperature,  =  -20,  -15,  -10°C; 

(ii)  Channel  slope,  =  8  x  10-5  and  2  x  10-^ 
(typical  for  pools  upstream  of  Hannibal); 

(iii)  Downstream  depth,  =  8, 1 2,  and  16  m  (nom¬ 
inally  26,  39,  and  53  ft,  respectively); 

(iv)  Pool  length,  £  =  26  to  162  km  (16  to  101  mi.) 

(v)  Composite  Manning’s  coefficient  n^.  =  0.033 
(based  on  representative  values  for  the  Ohio 
River  and  similar  large  rivers). 

TTie  numerical  results  were  obtained  for  certain 
combinations  of  Sq  and  Yf,  which  yielded  upstream  flow 
depths  greater  than  about  3  m. 

The  curves  presented  as  Q„  vs  t^  and  0^  vs  F  in  Fig¬ 
ures  19and  20,  respectively,  are  similarto  those  present¬ 
ed  in  Figures  15-18  for  the  Hannibal  and  Mon'mmery 
pods.  The  Q„  vs  V  curves  are  presented  for  =  '',0°C 
but,  as  explained  above  under  Ice  Formation  in  Hanni- 
baland  Montgomery  Pools,  are  approximately  valid  for 
Tg  =  -15  and  -10°C,  because  air  temperature  has  a 
relatively  insignificant  influence  on  V.  Variations  of  0^ 
with  Yf  for  different  pool  lengths  are  shown  in  Figures 
21  and  22,  with  Jp  =  8  X  10-5  and  2  x  IfH,  respectively. 
Values  of  0„  for  the  various  flow  conditions  were 
obtained  from  Figure  1 9.  Because  the  stabilit'’  criterion 
for  ice-cover  progression  (F  <  0.0  5)  does  not  depend  on 
air  temperature,  neither  do  values  of  0^.  Consequently, 
Figures  21  and  22  are  valid  for  all  air  temperatures. 

Air  temperature  may  in  fact  affect  0^,  but  in  ways  not 
accounted  for  in  the  present  simulation.  For  example, 
by  influencing  the  rate  and  quantity  of  frazil  ice  growth, 
it  might  affect  the  genesis  of  frazi  1  from  crystals  to  pans, 
so  that  all  frazil  formed  does  not  result  in  frazil  pans,  or 
pan  dimensions  may  differ  from  the  102-mm  ',4-in.) 
thickness  assumed  herein.  In  addition,  as  air  tempera¬ 
ture  influences  the  solidification,  or  subsequent  freez¬ 
ing,  of  an  accumulation  cover  as  it  forms,  it  could  affect 
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Figure  19.  Time  for  ice-cover  formation  for  S,,  =  ^  x  =  -20°C. 
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Figure  19  (cont'd).  Time  for  ice-cover  formation  for  =  2  x  /(H,  Ty  =  -20°C. 
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Figure  19  (cont'd).  Time  for  ice-cover  formation  for  S„  =  8x  10-^,  =  -20°C. 
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Figure  19  (cont'd).  Time  for  ice-cover  formation  for  S„  =  2  x  /CH,  Tg  =  -I5°C. 
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Figure  19  (cant’d).  Time  for  ice-cover  formation  for  =  5x  /O--^  T„  =  ~I0°C. 
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Figure  19  (cont'd).  Time  for  ice-cover  formation  forS„  =  2  x  10^,  =  -JO°C. 
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Figure  20.  Volume  of  ice  in  ice  cover  for  10-^,  =  -20°C. 
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Figure  20  (cont  ’d).  Volume  of  ice  in  ice  cover  for  =  2  x  10^,  Tj,  =  -20  °C. 
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cover  stability.  The  influences  on  of  air  temperature, 
however,  remain  a  matter  of  speculation  as  they  are  not 
yet  generally  documented. 

Values  of  increase  with  increasing  X/,  decreasing 
and  decreasing  t,  as  one  would  expect.  In  addition, 
Qo  is  very  sensitive  to  downstream  depths;  for  example. 
Qq  for  a  poo!  length  of  50  km  and  channel  slope  of  2  x 
10^  increases  from  150  to  520  mVs  (5,300  to  18,360 
ft^/s),  for  a  change  in  //  from  12  to  14  m.  If,  during  a 
period  of  ice-cover  formation,  it  is  permissible  to  in¬ 
crease  Y^,  significantly  higher  flows  could  be  used  to 
form  frazil-generated  ice  covers  for  reduced  periods, 
Also  keep  in  mind  that  lesser  values  of  would  be 
needed  for  ice  covers  to  form  over  portions  of,  rather 
than  entire,  pool  lengths. 


The  present  study  indicates  that  flow  regulation 
(primarily  discharge  regulation)  is  a  feasible  approach 
to  controlling  ice-cover  formation  on  pools  of  the  upper 
Ohio  River. 

CONCLUSIONS  AND  RECOMMENDATIONS 

To  minimize  frazil  ice  generation  and  thereby  pre¬ 
vent  the  occurrence  of  deleterious  frazil-ice-generated 
jams  in  rivers,  rapid  ice-cover  formation  is  required. 
One  way  to  achieve  this  is  by  means  of  flow  regulation, 
especially  discharge  regulation.  A  numerical  model  is 
developed  to  simulate  ice-cover  formation  by  frazil-ice 
accumulation.  The  simulation  model  is  applied  to  deter¬ 
mine  the  time  required  to  form  accumulation  ice  covers 
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Figure  21.  Optimal  river  discharge  for  =  5  x  10-^. 
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Figure  22.  Optimal  river  discharge  for  5^-2  x  1(H. 


28 


over  individual  river  pools  that  are  stage-regulated.  The 
numerical  results  of  the  simulation  model  indicate  that, 
for  individual  pools  or  reaches,  there  exists  an  optimal 
flow  discharge  that  corresponds  to  a  minimum  time  re¬ 
quired  for  the  ice  cover  to  form.  Optimal  discharge  in¬ 
creases  with  increasing  downstream  depth  of  flow,  de¬ 
creasing  pool/reach  length,  and  decreasing  channel  slope. 
In  addition,  time  to  form  an  ice  cover  decreases  with 
decreasing  air  temperature.  Figures  19  through  22  give 
information  on  times  and  ice  volumes  associated  with 
ice-cover  formation  over  several  sets  of  pool/reach  bed 
slope,  length,  and  downstream  depth. 

Based  on  an  example  involving  two  pools  of  the  Ohio 
River,  at  Hannibal  and  Montgomery,  it  appears  feasible 
to  minimize  the  amount  of  frazil  ice  generation  using 
existing  navigation  dams,  w!iich  would  serve  as  ice 
retention  structures  and  flow-stage  regulators,  and  exist¬ 
ing  flood-control  reservoirs  on  tributaries  that  could  be 
used  to  regulate  flow  discharge.  The  Ohio  River  and  its 
principal  tributaries,  the  Allegheny  and  the  Monongahe- 
la  rivers,  however,  comprise  a  fairly  complex  system  of 
several  pools,  not  just  a  single  pool.  Therefore  it  is  of 
interest  to  minimize  frazil  ice  growth  and  time  for  cover 
formation  over  extensive  portions  of  the  Ohio  River;  at 
least,  over  several  pools,  especially  in  the  vicinity  of 
potential  problem  locations.  In  this  regard,  the  problem 
of  controlled  ice-cover  formation  is  more  complex  for 
the  Ohio  River  than  it  is  for  the  St.  Lawrence  and  Lule 
River  examples  cited  in  Section  2,  Control  of  River  Ice 
Formation.  To  minimize  frazil  ice  growth  in  large  por¬ 
tions  of  the  Ohio  River,  it  is  necessary  to  treat  the  Ohio 
as  a  flow  system  such  as  indicated  in  Figure  6.  It  is 
beyond  the  scope  of  this  study  to  provide  a  detailed  guide 
as  to  how  flow  should  be  modified  or  regulated  along  the 
entire  Ohio  River.  That  task  is  a  logical  next  step  in  fully 
determining  the  feasibility  of  flow  regulation  to  control 
ice-cover  formation  on  the  Ohio.  It  would  require  a 
detailed  flow  model  that  simulates  a  series  of  pools, 
major  tributaries,  and  linked  flood-control  reservoirs. 

An  additional  logical  next  step  would  be  to  conduct  a 
field  experiment  aimed  at  verifying  the  proposed  strate¬ 
gy  of  flow  regulation  for  ice  management.  The  experi¬ 
ment  should  focus  on  a  single  pool.  A  reach  on  the 
Allegheny  River  below  Kinzua  Dam,  where  discharge 
could  be  controlled  closely  and  shore  ice  as  well  as  ice- 
cover  formation  could  be  monitored  in  detail,  is  recom¬ 
mended  as  a  possible  site  for  the  field  experiment.  In¬ 
formation  gained  from  such  a  field  experiment  would  be 
useful  in  further  development  of  a  simulation  model.  In 
particular,  it  would  provide  substantial  information  on 
the  genesis  of  frazil  ice  and  stability  criteria  for  ice- 
cover  progression,  important  two-  and  three-dimension¬ 
al  effects  not  accounted  for  in  the  present  one-dimen¬ 
sional  simulation. 


Further,  in  forming  accumulation  ice  covers  over 
critical  pools,  consideration  should  be  given  to  passing 
to  downstream  pools  ice  formed  on  upstream  pools.  In 
this  manner,  periods  required  for  ice-cover  formation 
might  be  reduced.  This  consideration  touches  on  an 
issue  frequently  confronting  operators  of  lock  and  dam 
installations:  whether  or  not  to  pass  broken  ice  con¬ 
veyed  by  the  river.  This  issue  has  to  be  addressed  in  the 
context  of  coordinated  efforts  to  either  form  or  release 
ice  covers  so  that  damaging  ice  jams  do  not  form. 
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APPENDIX  A:  EXPRESSIONS  FOR  COEFFICIENTS 

The  coefficients  a^,  in  eq  10  and  a,-,  bf,...,  in  eq  1 1  are  given  below: 

=  bi  =  -S~ 


=  =  andei  =  -5-(Gi+i-Gi) 

2at  Ax  Ax 

a,'= -^(Gw-Gi)(4) 

Ax  U2/i+, 


4Ax  L'A  /,+]  VA  /j  J 1 14  /j+j 


+  r;~[(3'i+l  -  >1)  fij+l  +  (Aj+I  +  Aj)] 
2Ax 


.  f  .  (1  -  p)s,„]  r)  j 


2Ax 

*f=  +  lQ\  +  Qi+i-Qi 

2At  AxlvA/i+i  VA/i  Aj^l 


2  '  i, 

c,= 

-  [{>’1+1  -  3'i)  -  Ui+I  +  Aj)] 

2Ax 

-  f  [[P^n  +  (1  -  (Ai^,  +  ^') } 


(i1i+i  -  Oi)  Bi 


d{=  -■J-+SS.  1^]  +  (Q|  -  gj+l-gj 
2Af  AxLW/i+i  'a/j  Aj 
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-“7“  Ui+i  +  ^^1)  (yi+1  -^i) 
2Ax 

2Ajc 

where  A:  =  iA/?2/3  andr  =&. 

«  dy 
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APPENDIX  B:  LISTING  OF  COMPUTER  PROGRAM 
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SIMULATION  01  THE  ICE  GENERATION  IN  AM  UMSTEASY  OPEN  CHANNEL 
SYSTEM. 

This  program  can  handle  a  ayatem  of  open  channels  with 
weirs,  storage  ponds. 

Currently  this  program  can  handle  only  one  node  point  with  several 
branched  channels  but  it  can  be  easily  «K>dified  to  handle  a  system 
with  several  node  points. 

Looped  channels  are  not  considered  in  this  program. 

All  channel  informations  should  be  given  in  order  of 
computational  sequence,  i.  e. , 

O/S  of  channel  1. - >  Node  point 

D/S  of  channel  2.  — — — — >  Node  point 

- > 

0/S  of  channel  n. - >  Node  point 

Rode  point  ______ - >  d/s  of  main  channel 


COMMON  /COl/ 

* 

COMMON  /C02/ 
COMMON  /C03/ 
COMMON  /C04/ 
COMMON  /C05/ 
COMMON  /C06/ 
COMMON  /C07  ' 
COMMON  /C08/ 
COMMON  /C09/ 
COMMON  /CIO/ 
COMMON  /Cll/ 
COMMON  /C13/ 
COMMON  /C14/ 
COMMON  /CIS/ 


A(34) ,B(34) .BW(34) .H(34) ,E{34) ,S(34) , 
CK(34),MCY(:'4),CKST(34/  Z(34) 
Q(34),V(34),Y(34) 

E(34),F(34),ITYP(34) 

YTDB(30J  ,QTDB(30)  ,QDB(30)  ,YDB(30)  ,NQDB,NYDB 
QT0B( 10,2) ,QDB( 10,2) ,NQ0B(2) ,HCHAN 
N,NT,TIME,I,KT,MNT1 ,LE 
DELY 

S0,ZL,ZR,1X>Y 

GRAY, SWI.DELTAT, ALPHA, BETA.THETA 
NN0D,N0D(2)  ,H0DU2) 

WY.WB.CKO 

DELX04),X(34) 

SORF.YBASE 

ICEFLG(34) ,T(34) ,C(34) ,aH(34) 


Following  are  the  key  variables  of  this  program. 

A(I);  area  ot  the  cross  section. 

B(I):  width  of  the  cross  section  at  the  top. 

BW(I):  width  of  the  cross  section  at  the  bottom. 

CKST(I):  Strickler's  cojfficient. 

CMU:  weir  coefficient. 

DELY;  dy. 

DPY:  dp/dy. 

H(I):  water  depth. 

P:  wetted  perimeter. 

R(l);  hydraulic  radius. 

S(I):  energy  slope. 

Q{I):  water  discharge. 

V(I):  velocity. 

Y(I);  water  surface  elevation  from  the  datum. 

B(l);  coefficient  E, 

F(I):  coefficient  F. 

YTDB(K):  time  at  which  the  D/S  water  surface  level  is  given. 
QTDB(K):  time  at  which  the  D/S  water  discharge  is  given. 
YDB(K);  given  D/S  water  surface  level. 

QTB(K):  given  D/S  water  discharge. 

NQDB:  number  of  diacharse  data  for  the  D/S  boundary  co_dition. 
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C  NYDB:  number  of  stage  data  for  the  0/S  boundary  condition. 

C  QTUB(K):  time  at  which  the  0/S  discharge  is  given. 

C  Ql)B(K):  given  U/S  discharge. 

C  NQ0B:  number  of  discharge  hydrograph  data  for  the  0/S  B. C. 

C  HNT:  number  of  time  steps  for  computation. 

C  NST:  number  of  time  steps  for  stabilisation. 

C  ICEFLG:  flag  for  ice  cover,  no  ice  “  0,  ice  covered  «  1. 

C  H:  number  of  space  intervals. 

C  HCEAN:  number  of  channels  branched  from  a  node  point. 

C  NNOD;  same  as  HCHAN. 

C  NOD(K):  point  No.  of  junction  point  of  Ktb  channel. 

C  NODKK);  point  No.  of  node  point  connected  to  current  junction 

C  point. 

C  NT;  time  index. 

C  TIME:  elapsed  time  from  the  starting. 

C  SO:  bed  slope. 

C  ZL:  slope  of  the  left  bank. 

C  ZR:  slope  of  the  right  bank. 

C  WY:  level  of  a  weir. 

C  VB:  width  of  a  weir. 

C  SURF;  surface  area  of  a  storage  pond. 

C  TBASE:  base  level  of  a  storage  pond. 

C 

C  ITYP-1  :  1-D  point 

C  2  :  D/S  point,  Q"Q(t) 

C  3  :  storage  pond 

C  4  :  junction  point 

C  5  :  weir 

C  6  ;  D/S  point,  Q^(y) 

C  ■/  :  D/S  point,  Q“Q<t) 

C  8  :  D/S  point,  y“y(t) 

C 

C 

DIMENSION  TYPECS) .PHAME(34) ,NAHE(34) ,XO(34) 

DIMENSION  HYDRO(28,200) ,FRONT(200) ,VOtICE(200J 

DIMENSION  IQHYD(  10)  .ITHYDf  10)  ,ITHYD(  10)  ,ICHYD(  10)  ,  IZHYD(IO) 

DIMENSION  TPT(30),ZHAME(1),KNAME(2) 

NCHAN=0 

NQDB=0 

NYDB'O 

NNOD=0 

READCl,*)  N.NOPT 

IF(N0PT  .EQ.  0)  GO  TO  810 

READd,*)  GQ.GY.GS.GDELX.GBH.GCKST 

X(1)-0.0 

Z(N)-0.0 

Y(N)-GY 

Ee-GDELX*1000.*GS 

H0-( GQ / ( GCKST*GBW*GS**0.5 ) )**0. 6 

DO  800  I-1,H 

ITYP(I)-1 

DELX(I)-GDELX 

BW(1)-GBW 

CKST(I)-GCKST 

Q(I)-<3Q 

IF(1  N)  GO  TO  800 

X(I+1)-X(I)-»DELX(I) 

Z(N-I)-ZiN-I-»l)-»Ee 
Y(N-I)-Z(H-I)+H0 
1F(Y(N-I)  .LT.  GY)  Y(H-I)-GY 
800  CONTINUE 
810  DO  10  I-1,N 

READCl,*)  PNAME(I),ITYP{I),X(1),M;lX(I),Z(I),BW(I),CKST(I) 
ZNAME(1)-PNAME(I) 

DECODEC  1, 700, ZNAKE)  HAMECD 

WRITEC6,610)  I,PNAMECI),ITYPCI),ICI).DEI.XCI),ZCI).BWCI),CKSTCI) 

DELXCI)-DELXCI)*1000.0 

XOCD-XCI) 
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X(I)-X(I)*1000.0 

GO  TO  (11,2,3,4,5,6,6,6),  ITTPd) 

C 

C  D/s  of  a  channel  :  input  dlacharge  hjrdrograph  it  required  aa 
C  the  boundary  condition. 

C 

2  NCHAN-NCBAN^l 
KNAME (NCH4H ) -NAHE ( I ) 

HEAD( 1 ,*)  HQ, (QTDB( J,NCHAH) ,QDB( J.HCEAH) , J-1 ,»Q) 

NQDB(HCHAH)-Hq 

WRITE(6,615)  (QTDB(J,HCaAH),QDB(J,HCHAH),J-l,»Q) 

GO  TO  11 
C 

C  Storage  pond  :  surface  area  in  hectare  and  base  level  are 
C  required. 

C 

3  READd,*)  SDRF.YBASE 
WRITE(6,620)  SDRF.YBASE 
SDRF-SDRf*10000.0 

GO  TO  10 

4  NN0D>’NN00-f.' 

N0D(NN0D}-I 
READd,*)  NODUmiOD) 

WRITE(6,630)  NODUHNOD) 

GO  TO  10 

C 

C  Veir  :  weir  elevation,  width,  and  weir  coefficient  are  required. 

c 

5  READd,*)  WY,WB,CM0 
HRITE(6,640)  HY,»B,aro 
GO  TO  10 

6  IF<ITYF(I)-7)  7.8,9 
C 

C  D/S  B.  C.  ;  Rating  curve  :  Q"Q(y) 

C 

7  WRITE (6, 650) 

GO  TO  10 

C 

C  D/S  B.C.  ;  Q-Q(t) 

C 

8  READd,*)  NQDB,(qTDB(J),QDB(J).J-l,HQDB) 

WRITE(6,660)  (QTDB( J) ,QI»( J) , J-1,HQDB) 

GO  TO  10 

C 

C  D/S  B. C.  ;  y“y{t) 

C 

9  READCl,*)  MYDB,(YTDB(J),YDB(J),J-1,HYDB) 

WRITE(6,670)  (YTDB( J) ,YM( J) , J-l.HYM) 

GO  TO  10 

11  CONTINDE 
10  COHTIHDE 
C 

INITIALIZE  THE  VARIABLES  AND  READ  INITIAL  CONDITIONS 

SO»(Z(N-l)-Z(N))/DELX(N-l) 

IF(NOPT  .EQ.  0)  RKADCl,*)  (Y(I) .Qd) ,I-1,N) 

READd,*)  GRAY, ALPHA, BETA, THETA 
READd,*)  TYPE,ZL,ZR 
READd,*)  THAX,DELT,NST 
DELXAT-DELT*60.0 
MNT-XHAX*60. 0/ DELTAT+0. 01 
MNTl-MNT+l 
NTT-MNT+NST 

READ  TIMES  FOR  SIRFACE  PROFILE. 

READd,*)  NIPT,(TPT(J),J-1,NIPX) 
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READ  SECTION  NO.  FOR  HYDROCRAPH. 

READd,*)  NQHYD,(IQHyD{I),I-ltNQHYD) 

READd,*)  NYHYD.dYHYDCD.I-l.NlfHYD) 

READCl,*)  NTHYD,(ITHYDd),I-l.HTHYD) 

READd,*)  NCHYD,dCHYD(I),I-l.HCHYD) 

REaD(  i  ,*)  NiuiD.dZHYDd)  .I-I.MZHYD) 

NQL-NQHYD+1 

NYF-NQL+l 

NYL-NQL+HYHYD 

NTF-NYL+I 

NTL-NYL+NTHYD 

NCF-NTL+1 

NCL-NTL+NCHYD 

NZF-NCL+1 

NZL-NCL'fNZBYD 

PRINT  INITIAL  VALDES 

WRITE(6.200)  N.NCHAN,<atAV.TYFE,SO.ZL.ZR,TMAX.IXELTAT, 

#  MNT.NST.ALPHA, BETA, THETA 

BEGIN  LOOP  FOR  TIME 

TIKE-0.0 

NT-1 

CALL  THICON 
25  DO  30  I-1,N 
30  Hd)-Yd)-Zd) 

KT-NT-NST+1 
CALL  ICE 
CALL  SECTN 
DO  40  I-1,M 

40  V(I)-Qd)/A(l) 

IF(HT  .LT.  HST)  CO  TO  46 

COMPUTE  TEMPERATURE  DISTRIBUTION  AND  ICE  CONCENTRATION. 

CALL  TEMP 
CALL  ICECVR 
CALL  CONDUC 

I>(TH(N)  .GT.  0.01  GO  TO  35 
FR0NT(KT)-0.0 
V0LICE(KT)-0.0 
GO  TO  38 

35  SDMICE-0.0 
Nl-N-1 

DO  36  I-1,N1 
NI-N-I 

IF(TH(NI)  .LE.  0.0)  GO  TO  37 
8UMICE-SUMICE+TH(NI)*DELX(NI)*B(BI) 

36  CONTINUE 

37  VOLICE(KT)-SOMICE 
IRONT<KT)-(X(N)-X(NI+1))/1000. 

CREATE  OUTPUT  FILE 

38  HYDR0d,Kr)-TIME/60,0 

DO  41  J-l.NQHYD 

41  HYDRO(J+l,Kr)-qdQHYD(J)) 

DO  42  J-1,NYHYD 

42  HY1»0(J+NQL,KT)-Y(IYHYD(J)) 

DO  43  J-1,HTHYD 

43  HYDRO(J+NYL,KT)-T(ITHYD(J)) 

DO  44  J-1,NCHYD 
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44  HYDR0(J+IJTL,KT)»C{1CHYD(J)) 

DO  45  J-l.HZHYD 

45  HYDRO(J+NCL,KT)-TH(lZHYD(J)) 

46  DO  49  J-1,NIPT 

1F(ABS(HKE-WT(J>)  .GT.  l.E-2  .OE.  KT  .IT.  NST)  GO  TO  49 
miITE(6,300)  TIKE 
DO  48  K-1,NCHAN 
WRITE (6, 3 10)  K 
DO  47  I-1,H 

IF(NAMB(I)  .ME.  KNAME(K))  GO  TO  47 

W&ITE(6,320)  I.PNAME{I).XO(I),Y(I).H(l),A(I),Q(I),T(I),C(I),ra(I) 

47  C(»ITIMDE 

48  OIHTIMOB 

49  COMTIMDE 

IF(NT  .GE.  MTT)  GO  TO  100 
IF(NT  .LT.  MST)  GO  TO  50 
TIME-TIME +IKLT 

50  CALL  DSWEEP 
HT-NT+1 

GO  TO  25 

100  MRITE<6,500)  tPNAME(IQHYD(  J)  )  ,J-1  ,BQHYD) 

DO  110  J-1,MMT1 

110  WRITE(6,510)  HYDR0(I,J),(HYDR0(I,J).I-2,1IQL) 

WRITE(6,520)  (PNAME(IYHYD(J)),J-1,HYHYD) 

DO  120  J-1,MNT1 

120  WRITE(6,510)  HYI»0(1.  J)  ,(HYDRO(I,  J)  ,I-HYF,MYL) 

WRirE(6,530)  (PNAMEdTHYDC J)) , J-1,KTHYD) 

DO  130  J-l,MNTl 

130  WRITE(6.512)  HYDRO(l,  J)  ,(BYDRO(I,  J)  ,I-MTF,MTL) 

WR1TE(6,540)  (PMAME(ICHYD(J))  ,J-1,MCHYD) 

DO  140  J-1,MNT1 

140  WRITE(6,5H)  HYCR0(1,  J)  ,(HYDRO(I,  J)  ,I-HCF,HCL) 

WRITE(6,550)  (PMAME(IZHYD(J)),J-1,1IZHYD) 

DO  150  J-1,MNT1 

150  WRITE(6,512)  HYDRO(l,J)  ,(HYI»lO<I,  J)  ,1-NZF,HZL) 

WRITE(6,560) 

WRITE(6,5i3)  1HYDR0(1,J),FR0NT(J),V0LICE(J),J-1,MNT1) 

STOP 

200  F0RMAT{1H1.//5X, '*♦****♦  PROBLEM  DESCRIPTIOM** 
,/// ,'  <  GEOMETRICAL  CONDITIOH  >',//, 

Number  of  points  -  ',13,/ 

Ntimber  of  channels  -  ',13,/ 

Gravitational  acceleration  »  ',F5.2,'’(m/sec**2)',/ 

#'  Channel  type  -  ',8A4,/ 

#'  Bed  slope  at  D/S  -  ',F8,6,/ 

#'  Slope  of  left  bank  -  1  ;  ',F3.1,/ 

#'  Slope  of  right  bank  -  1  :  ',F3.i,/// 

#'  <  NDMERICAL  CONDITION  >',// 

Maximum  time  -  ',F7.1,'(sd.n)',/ 

#'  Time  step  -  ' ,F7.1,'(sec)',/ 

#'  No.  of  time  step  for  unsteady  flow  -  ',13,/ 

#'  No.  of  time  step  for  stabilization  -  ',13,/ 

#'  Energy  correction  coefficient  -  ',F3.1,/ 
f'  Weighting  coefficient  in  space  •  ',F4.2,/ 

#'  Weighting  coefficient  in  time  -  ',F4.2,//) 

300  FORMAT (IHl,//,'  Time  -  ',F7.1,'(min)') 

310  FORMAT (//,1 OX, 'Channel  No.  -  ',12,// 

#  2X,'No.  Name  Z(km)  Y(n)  h(m)  A(m'*te}  Q(cms)' 

TEMP(C)  C(X>  ICE(m)',//) 

320  FORMAT  C 14 ,3X, A4 ,  F6. 1 ,  F8. 3 ,  F7.  3 ,  F8. 2  ,F8. 2 ,  F9.  4 , 2PF7. 3 ,0PF8. 4) 

500  F0RMAT(1H1.//,5X,'<  DISCHARGE  HTDR0CRAPH>' 

#  ,//,14X, 'Water  Discharge  (cms)',// 

#  ,3X.'TIME'.4X,10U4,5I),//) 

510  F0RMAT(lX,r6.2,10F9.3) 

511  FORMAT<lX,0PP6.2,2P10F9.3) 

512  FORMAT(1X,F6.2,10F9.4) 

513  FORMAT(1X.F6.2,F10.2.F12.0) 
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520  FORMAT(lHl,//.5X,'<  STAGE  HYDROGRAPH  >',// 

#  tl4X, 'Water  Surface  Level  (m)',// 

#  ,3X,'TIME',4Z,10(A4,5X),//) 

530  FORMAT  (IHl.//, 5  X,'<  TEMPERATDRE  DISTR  >',// 

#  t 14X, 'Average  temperature  (C)',// 

#  ,3X,'TIME'.4X,10(A4,5X>,//) 

540  F0RMAT(1H1,//,5X,'<  ICE  COPCENTRATIOH  >',// 

#  il4X, 'Average  concentration  (.X)',// 

#  ,3X,'TIME'.4X.10{A4.5X),//) 

550  F0RMAT(1H1,//,5X,'<  ICE  THICKNESS  >',// 

#  ,14X,'Ice  cover  thickneas  (m)',// 

#  ,3X.'TIME',4X,10U4,5X),//) 

560  FORMAT(lHl,//,5X,'<  ICE  ACCOMDLATION  >',// 

#  .3X,'TIME',4X,'X  (kn)',3X,'V0L  (m**3r,/) 

610  FORMAT<//,1X,60( IE. 'Ho.  ',I2,2X,A4,2X,'Type-' ,I2,2X 

#  ,'I-',P4.1,2X,'])X-',F4.1,2X,'Z-',F5.2,2X,'BW-',F5.1 

#  ,2X,'Kat--,F4.1./,lX.60('--)) 

615  FORMAT(/,5X,'<  OPSTRKAM  BODNDAEY  POINT  >',// 

#  ,5X,'***  Input  Hydrograph  ***',// 

#  T(min)  Q(cm8)',/,(F11.3,5X,F11.3)) 

620  F0RMAT(/,5X,'<  STORAGE  POND  >',// 

#  .SX, 'Surface  Area  -  ',F4.1,'(ha)',/ 

#  ,5X,'Bottom  Elevation  «  ',F5.2,'(ia)') 

630  FORMAT(/,5X,'<  JUNCTION  POINT  >',// 

#  .5X,'  D/S  point  -  No. ',12) 

640  P0RMAT(/,5X,'<  RECTANGDLAR  WEIR  >',// 

t  ,5X, 'Elevation  of  the  Weir  «  ',F4.1 ,'(m) ',/ 

#  ,5 X, 'Width  of  the  Weir  ■=  ',F4.1,'(m)',/ 

t  ,5X,'Weir  Coefficient  •  ',F4. 1) 

650  P0RMAT(/,5X,'<  D/S  BODNDARY  POINT  >',// 

#  ,5X,'***  Rating  Curve  ***-) 

660  FORMAT(/.5X,'<  D/S  BOONDiRY  POINT  >',// 

#  ,5X,'***  Discharge  Hydrograph  ***',// 

#  T(nin)  Q(cm8)',/,(F11.3,5X.F11.3)) 

670  F0RM&T(/.5X,'<  D/S  BOOHaiRY  POINT  >',// 

#  ,5X,'***  Stage  Bydrograph  ***',// 

#  T(iBin)  Y(m)',/,(F11.3,5X,F11.3)) 

700  FORMAT  (AU 

END 


SUBROUTINE  DSWEEP 

COMMON  /COl/  A(34),B(34),BW(34),H(34),R(34).S(34), 
*  CK(34),DKY(34),CXST(34),Z(34) 

COMMON  /C02/  Q(34),V(34),T(34) 

COMMON  /C03/  E(34),F<34),ITYP(34) 

COM»K)N  /C06/  H,NT,TIME,I,Kr,MNTl,LE 
COMIBIN  /C07/  DELY 

COMMON  /C09/  GRAY, SWI, mXAT, ALPHA, BETA.THEXA 
COMMON  /CIO/  NNOD,HOD(2),NODU2) 

COMMON  /C12/  AA,BB,CC,DD.GC,AA1,BB1,CC1,0D1.GG1 
COMMON  /C13/  DELX(34),X(34) 

COMMON  /C15/  ICKFLG(34).T(34).C(34),TH(34) 
DIMENSION  CL(34),CM(34).CN(34) 

IF(NT  .GT.  1)  GO  TO  5 

G2-GRAV*0.5 

GT-GBAV*THETA 

GT2-GT*0.5 

R2T-0.5/DELTAT 

H1..N-1 

5  DO  100  1«1,N1 
DELTAX-DELX(I) 

GO  TO  (10,20,30,40,50),  ITYP(l) 

20  CALL  UPBC 
GO  TO  60 
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30  CALL  FOND 
GO  TO  90 
40  CALL  JONC 
GO  TO  100 
50  CALL  WEIR 
GO  TO  90 

10  IF(DELTAX  .EQ.  DELX(I-ll)  GO  TO  70 
60  ADX>ALPHA/DELTAX 
A4DX»ADX*0. 25 
ATDX=AI)X*THETA 
AT2DX-ATDX*0.5 
AT4DX=AT2DX*0.5 
GT2DX-GT2/DELTAX 
G2DX-<;2/DELTAX 
70  QMQ-Q(I+1)-<J(I) 

VPV-V(I+1)+V(I) 

AMA-A(I+1)-A(I) 

TMT-T(I+1)-Y(I) 

TMT-(TH(I+1)-TH(I))*SWI 

AFA-A(I+1J+A(I) 

SPS-BETA*S(I)+(1.-BETA)*S(I+1> 

QA1-V(I+1)/A(I+17 

QA-V(I}/A(I) 

C 

C  COEFFICIENTS  OF  THE  DYNAMIC  EQUATION 
C 

AA«-ATDX*QMQ*QA1*B(  I+l)  -AT4DX*VPV*(  VPV-2.  *AMA*<QA1 }  i*B  ( I+l ) 

#  +GT2DX*({YMY+TMr)*B(I+l)+APA> 

#  •tGT2*(SPS*B(I+l>-2.*APA*(l.-BETA)*S(I+l)/CK(I+l)*DKY(I+l)) 
CC-ATDX*QM!}4QA*B(I)-AT4DX*VPV*(VPV+2.  *AMA«QA)*B(  I) 

#  -GT2DX*(<YMY+TKr)*B(I)-APA) 

#  -GT2*(  SPS*B(  I)  -2,  *APA*BETA*S(I  l/CKd  )*DKY  ( I )  ) 

BB-R2T+ATDX*  (  VPV+QMQ/A(  I+U  )  -AT2DX*VPV*AMA/A  ( I+l ) 

t  +GT*APA*(l.-BETA)*ABS(Q(Id))/{CK(I+l)*CK(I+l)) 

DD— R2T+ATDX*(VPV-QMQ/A(I»+AT2DX*VPV*AMA/A{I) 

#  -GT*APA*BETA*ABS(q(I))/(CK<I)*CX(I)) 

GG— ADX*VPV*QMQ+AADX’*VPV'»VPV*AMA-G2DX*APA*(YMY+TKr)-G2*APA*SPS 
C 

C  COEFFICIENTS  OF  THE  CONTINUITY  EQUATION 
C 

AA1-R2T*B(I+1) 

CCl— R2T*B(1) 

BBl-THETA/DELTAX 

DDl-BBl 

GGl— QMQ/DELTAX 
C 

90  CDCD-CC1*DD-CC*DD1 

CL( I ) -( AA1*DD-AA*DD1) /CDCD 
CM(  I )  -(  BBl  *DI>-BB*DD1 )  /CDCD 
CN ( I ) -( DD 1*GG-DD*GG1 ) /CDCD 
CDE-CC1+DD1*E(I) 

BCC-BB1-CM(I)*CDE 
E(I-fl ) -(  CL(  D^CK-AAl ) /BCC 
F(I+1)-(CN(I)*CDE>DD1*F(I)+GG1)/BCC 
100  CONTINUE 
CALL  DNBC 

DELQ=E(N)*DELY+F(N) 

Y(N)->Y(N)+DELY 
Q(N)-Q(N)+DELQ 
DO  200  J-1,N1 
I-N-J 

IFdTYPd)  .NE.  4)  GO  TO  160 
C 

C  JUNCTION  POINT 
C 

DO  150  K-l,NNOD 

IFd  .BQ.  NOD(K))  CO  TO  155 
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150  CONTINUE 

155  DELT-Y(NODl(K))-T(I) 

GO  TO  165 

160  DELY-a(I)*DELY+CM<I)*DELQ+CN(I) 
165  Y(l)-Y(I)+IffiLY 

DEIXl-E(I)*DEiY+F{I) 

200  q( !)•<)(  I 
RETURN 
END 


<iniioAnTTMP 

COMMON  /COl/  A(34),B(34).BV(34},H(34),R(34),S(34}. 
#  CK(34),Drr(34),CR8T(34),Z(34) 

COMMON  /C02/  Q<34),V(34),Y(34) 

COMMON  /C06/  N.MT.TIMB.l.Er.MNTl.LE 

COMMON  /C08/  SO,ZL,Zi,WY 

COMMON  /CIS/  ICEILG(34),T(34).C(34}.‘tH(34) 

IF<IIT  .GT.  1)  GO  TO  10 

ZRL-ZR>ZL 

SL-SQRT(1.0+2L*ZL) 

SR-SQRTCl.O-fZR^ZR) 

IffY-SR+SL 
10  DO  20  I'l.H 

B(I)-BW(I)+H(I)*ZBL 
A(I)-(BW(I)+B{I))*H(I)*0.5 
P-BW(I)+H{I)*DPY 
IFdCEiTGd)  .BQ.  1)  P-P+Bd) 

Rd)-Ad)/P 

CKR23-CKSTd)*Rd)**0.  6666667 

I*Yd)-0.6666667*CKR23*(2.5*Bd)-Rd)*DPY) 

CKd)-CKR23*Ad) 

20  Sd)-qd)*ABS(Qd))/(CKa)*CKd)) 

RETURN 

END 


SUBROUTINE  WBC 

COMMON  /COl/  A(34),B(34),BW(34),H(34),R(34),5(34), 

*  CK(34),DKY(34),CKST(34),Z(34) 

COMMON  /C02/  Q(34),V(34),Y<34) 

COMMON  /C03/  E(34),F(34),ITYF(34) 

COMMON  /C04/  YTDB(301,QTDB(30),QDB(30;,YI»(30),NQDB,HYia 
COMMON  /C06/  N,NT,TIME,I,KT.MNT1.I^ 

COMMON  /C07/  DELY 
COMMON  /C08/  SO.ZL.ZR.EPY 

IF(NT  .GT.  1)  GO  TO  5 
CKS5»CKST(H)*S0**0.5 
NQDBl-NQDB-l 
NYDBl-NYDB-1 

5  IPdTYP(N)-7J  30,10,20 

INTERPOLATE  THE  DISCHARGE  BYDROGRAPH 

10  IFdlME  .GE.  QTDBd})  GO  TO  12 
DELY-(qDB( lJ-q(N)-F(B))/E(H) 

RETURN 

12  DO  15  K-l,NqDBl 

1F(. MOT.  (TIME  .GE.  qTl»(K)  .AND. 

*  TIME  .LT.  qTIS(K+ll})  GO  TO  IS 
DELY-(qDB(K)+(TIME-qTI»(K))*(qDB<K+ll-qDB(K)) 

*  /(qTDB(K+lJ-qTI»(R))  -q(H)  -F(N))  /  K(N) 

RETURN 

15  CONTINUE 

DELY-(qDB(NqDB  ) -<l(  N)-F(H)  ) /E(B) 

RETURN 
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c 

C  INTERfOLATE  THE  STAGE  BTD&OGRAFH 
C 

20  IF(TI1!E  .GE.  TTMCl))  GO  TO  22 
DELY-TDB( l)-T(N) 

RETUEN 

22  DO  25  E-l.mrDBl 

IF(. NOT.  (TIME  .GE.  TTia(K)  .AND. 

*  TIME  .LT.  YTDB(K*1)))  GO  TO  25 
DELY-YDB(K)+(TlME-TTia(E))*(TDB(K+U-TDB(K)) 

*  /(YTDB(K+l)-YTia(K))  -Y(N) 

RETDEM 

25  COMTINOE 

DELY-YDB(NYDB)-Y(H) 

RETURN 

C 

C  FLOW  IS  LOCALIY  UNIFORM  :  MANNING'S  EQ 
C 

30  CKYS-CKS5*R(N)**0.6666667 

DELY-{-Q(N)+CKYS*A(N)-F(N))/{E(N)-(1.666667*B(N) 

*  -0.6666667*R(R)*DPY)*CKYS) 

RETURN 

END 


SUBROUTINE  UFBC 

COMMON  /C02/  Q(34),V(34),Y(34) 

COMMON  /C03/  E(34).F(34),ITYP(34) 

COMMON  /C05/  QTUB(10,2),QDB(10,2),NQDB(2),NCHAN 
COMMON  /C06/  N,NT,TIME,I,ET,MNT1,LE 


DATA  NCH/1/ 

C 

C  INTERPOLATE  THE  DISCHARGE  HYDROCRAPB 
C 

1  NQDB1-NQDB<NCH)-1 
E(I)-0.0 

IFCTIME  .GE.  QTUBd.NCH))  GO  TO  5 
F(I)-<10B(1,NCH)-Q(I) 

GO  TO  20 

5  DO  10  K-l,RqDBl 

IF(. HOT.  (TIME  .GE.  QTUB(K,NCH)  .AND. 

#  TIME  .LT.  qTUB(K4^1,NCB)})  GO  TO  10 
P(I)-<J0B(K,NCH)+(TIME-<JTUB(K,NCH))*(QDB(K+1,HCH)-QDB(K,NCH)) 

#  /(QTDB(K+l,HCH)-QTDB(E,NCa))  -<}(I) 

GO  TO  20 

10  CONTINUE 

F(I)-QDB(NQDB(NCB}  ,NCB}-<2(I) 

20  NCH-NCE-i-l 

IF(NCH  .6T.  NCBAN)  NCH-1 

RETURN 

QID 


SUBROUTINE  FOND 

COMMON  /C02/  <l(34),V(34),Y(34) 

COMMON  /C06/  N,NT,TIMB.1.EI,)«T1,LE 

COMMON  /C09/  GRAV.SHl.IKLTAX.ALPBA.BETA.THEXA 

COMMON  /C12/  AA.BB,CC.D0,GG,AA1,BB1,CC1,DD1,GG1 

COMMON  /C14/  8URF,YBASE 

AAl-0.0 

BBl— THETA 

IF(Y(I)  .LE.  YBASE)  GO  TO  5 

CCl-SURF/DELTAT 

GO  TO  10 
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5  CCl-0.0 
10  DDl—THETA 

GGl-Qd+D-Kia) 

AA>1.0 

BB=0.0 

CC-1.0 

DD-0.0 

GG-Y(I)-Y(I+1) 

RETURN 

END 


SUBROUTINE  WEIR 

COMMON  /C02/  Q(34),V(34),Y(34> 

COMMON  /C06/  N.NT.TIME.I.KX.MNTI.LE 

COMMON  /C09/  GRAY, SWI.DELTAT.ALPHA, BETA, THETA 

COMMON  /Cll/  Wf.WB.CMU 

COMMON  /C12/  AA,BB,CC,DD,GG,AA1,BB1,CC1,DDI.GG1 
IF(NT  .GT.  1)  GO  TO  5 
CMBG»CIflJ*WB*<  2.0*GRAV)**0.5 
CMBG3-CMBG*0.57735 
CMBGE>CMBG*10.0 
S  AAl-0.0 
BBl-1.0 
CCi-0.0 
DDl-1.0 

GG1-<}(I)-Q(I*1) 

BB-0.0 

DD-1.0 

1F(Q{I)  .LT.  0.0)  GO  TO  10 
C 

C  POSITIVE  Q 

C 

TUS»Y(I) 

YDS-Y(I+1) 

GO  TO  20 
C 

C  NEGATIVE  Q 

C 

10  YUS-T(I+1) 

TDS-Y(l) 

20  TUSW-YDS-WY 
YDSW-YUS-HY 
YUSDS-YUS-YDS 

IF(YDSW  .GT.  0.6666667nUSW)  GO  TO  30 
C 

C  FREE  FLOW 

C 

IjqOS-CMBG3*TOSW»*0.5 

DQOS-0.0 

QW-0. 6666667*DQ0SnU5W 
GO  TO  SO 
C 

C  FLOODED  WEIR 

C 

30  IF(ABS(YDSDS)  .LT.  0.01)  60  TO  40 
YDSDS5-YUSDS**0.5 
T-0.5*YDSW/YDSDS5 
DQUS-CMBGn 
DQDS-CMBG*(YDSDS5-T) 

QW“CMBG*YUSDS5  *YDSW 
GO  TO  so 

40  DQDS-CMBGE*YDSW 

IIQDS-CMBGE*(YDSDS-YDSW) 

QW-DQDS*YUSDS 

SO  IF{Q(I)  .LT.  0.0)  GO  TO  60 
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MxCQDS 
CC—DQOS 
GG»Q(I)-QW 
EETURN 
60  AA— DQOS 
CC-DQDS 
GG-<}(1)+QW 
B£XURN 
END 


SDBRODTIME  JUNC 

COMMON  /C02/  Q(3A),V(3A),Y(3A) 

COMMON  /C03/  E{3A).F(3A),ITTP(3A) 

COMMON  /C06/  N,NT,TIME,I,Kr,MNTl,LE 
COMMON  /CIO/  MN0D,N0D(2).M0D1(2) 

DATA  EJtINC,FJDNC/2*0.0/ 

5  DO  10  K»1,HN0D 

IF(I  .EQ.  HOD(K))  GO  TO  15 
10  CONTINOE 
15  EJONC-EJDNC+E(I) 

FJDNC-FjnNC+Q(I)+E(I)*(Y(NODl(K))-Y(I))+F(I) 

IF(R  .NE.  NNOD)  RETDKN 

E(N0D1(K))-EJIINC 

F(NODKK)  )  -FJDHC-Q(NODllK)  ) 

EJUNOO.O 

FJDMC«0.0 

RETURN 

END 


10 

20 


30 


SUBROUTINE  ICE 

COMMON  /C06/  N,NT,TIME,I,Kr,MNTl,LE 
COMMON  /CIS/  ICEnX3(3A),T(3A),C(3A),TH(3A) 
IFiST  .GT.  1)  GO  TO  20 
DO  10  I-1,N 
iCEnx;(i)-o 

IF(TH(1)  .GT.  0.0)  ICEJlGd)-! 

CONTINUE 

IF(KT  .LT.  1)  RETURN 
DO  30  I-1,N 

IF(TH(I)  .GT.  0.0)  ICEFLGd)-! 

CONTINUE 

RETURN 

END 


SUBROUTINE  TEMP 


VERTICALLY  AVERAGED  WATER  TEMPERATURE  AND  ICE  CONCENTRATION 
ARE  COMPUTED  USING  THE  CHARACTERISTIC  BASED  B-P  ATH  ORDER 
INTERPOLATION  SCHEME. 


COMMON  /COl/  A(3A),B(3A),BW(3A),H(3A),H(3A),S(3A), 

*  CK(3A),DKY<3A),CKST(3A),Z(3A) 

COMMON  /C02/  Q(3A),V(3A),Y(3A) 

COMMON  /C06/  H,NT,TIME,l,Kr,MNTl,LE 

COMMON  /C09/  GRAV, SHI. mTAT, ALPHA, BETA.THETA 

COMMON  /C13/  DELX(3A),X(3A) 

COMMON  /C15/  ICEFLG(3A),T(3A),C(3A),TH(3A) 

COMMON  /C16/  TA{200),TTDB(200),TO(3A).CTDB(200),CO(3A) 
COMMON  /C17/  ROU.HWA.CLAINT 

DIMENSION  HO(3A)  ,HXO(3A)  ,HX(3A)  ,V0(3A)  .VXOOA)  ,VX(3A) 
DIMENSION  ALP0(3A) ,ALP(3A) ,TX0(3A)  ,TX(3A) ,TT(200) 
DIMENSION  CXO(3A),a(3A),CT(200) 

IF(KT  .GT.  1)  GO  TO  120 
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K£AO  AND  DEFINE  THE  JNIIIAL  VAIDES. 

READd,*)  HHA.ROT.SPHT.CLATNT 
WRITE (6, 5 00)  IKA,ROD.S?HT,CLATNT 
HWA>BWA*10000.  /  (  24.  *3600.  ) 

ROU-ROD*l.E6 

CALL  AIRTHP 

CALL  TICOND 

CALL  TUPBC 

CALL  Cl  CON  D 

CALL  COPBC 

CALL  DIFEX(TO,TXO,N) 

CALL  DIFFT(TTDB,TT,MNT1> 

Nl-N-I 

DO  100  I-1,N1 

IFdCEFLGd)  .HE.  0)  GO  TO  105 
100  CONTINDE 
105  LE-I 

CALL  DIFFX(C0.CZ0,N) 

CALL  DIFFT(CTUB,Cr,MKTl) 

CST-HWA/(ROtJ*SPHT) 

EPS-BWA/(ROU*CLATNT} 

DTAT2-DELTAT*0.5 

EPST2-EPS*DTAT2 

110  DO  111  I-1,N 
HOd)-Hd) 

VOd)-Vd) 

Td)-T0d) 

TXd)-TX0d) 

cd)-cod) 

CX(I)-CXOd) 

111  ALPOd)-CST/HO(I) 

CALL  DIFFX(HO,HXO,N) 

CALL  DIFFX(VO,VXO,N) 

RETURN 

120  CALL  DIFFX(H.BX,N) 

CALL  D1FFX(V,VX,H) 
m-KT-l 
DO  122  I-1,N 
122  ALPd)-CST/H(I) 

T(1)-TTDB(KT) 

TX( 1 ) -( ALP( 1 ) *(TA(KT )-T( 1) )-TT(KT ) ) /V( 1) 

TX(1)— TT(KT)/V(1) 

C(1)-CTUB(RT) 

CX( 1) — (EPS*TA{XX)/H( 1)+Cr  (XT) ) / V( 1) 

NOICE'LE 

IFdCEFLG(LE)  .NE.  1}  NOlCE-LE-fl 
DO  200  1*2 .NOICE 

COMPUTE  THE  TRAJECTORT  OF  THE  CHARACTERISTIC  LINE  USING  THE 
TIME  AVERAGED  WATER  VELOCITIES. 

XS I-X  (I ) -DTAT2*  (  V(  I )  ♦VOd-1 )  ) 
mXSI  .LT.  0.0)  GO  TO  150 
ITR-1 

125  DO  130  L-2,I 

IF(XSI  .GE.  X(L-l)  .AND.  XSI  .LE.  X(L))  CO  TO  135 
130  CONTINUE 
WRITE(6,300) 

300  FORMAT (///,1X,20('*'),'XSI  IS  OUT  OF  BOUNDART') 

STOP  200 
135  Ll-L-1 

QCOH-DrAr2*  (  V0(L)  -VO(U  )  )  /DELX(  LI ) 

XSI-(Xd  )  -DTAT2*(  Vd)4VO(L)  )  •K)CON*X(L)  )  /(QCON+1 . 0) 

1F(XSI  .LT.  0.0)  GO  TO  150 

1F(XSI  .GE.  X(L1)  .AND.  XSI  .LE.  X(L)  .OR.  ITR  .GE.  3)  GO  TO  140 

ITR-ITR+1 

GO  TO  125 
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c 

C  TRAJECTORY  OF  THE  CHARACTERISTIC  LIKE  FALLS  WITHIN  THE  CHANNEL 
C  REACH. 

C  INTERFOLATIOH  ON  SPACE  IS  NEEDED. 

C 

140  TA0-(X(L)-XSI)/DEUC(L1) 

TAD2-TAD*TAD 

COMPUTE  THE  COEFFICIENTS  OF  THE  4TH  ORIER  H-P  INTERPOLATION 
SCHEME  BASED  ON  THE  BESMITE' QUBIC  POLYNOMIAL. 

A1»TAD2*(3.-2.*TAD) 

A2=1.-A1 

A3=TAU2*(  1.  -TAH)*DELX(  LI ) 

A4=-TAD*(  1.  -TAD)**2*DELX(U  ) 

B1-6.0*TAD*(TAD-1. )/DELX(Ll) 

B2— B1 

B3-XAD*(3.*TAU-2.) 

B4-(TAD-1. )*(3.*TA0-1. ) 

IF(TO(I)  .LE.  0.0)  GO  TO  145 

LINEAR  INTERPOLATION  FOR  THE  FOLLOWING  QUANTITIES. 

143  ALPXSX l.-TAD)*ALPO(L)+TAD*ALPO(Ll) 

VXSI«=(  1.  -TAU)  *VO(  L)  +TAD*VO(  LI ) 

VXXSI-(  l.-TAU)*VXO(L)+TAU*VXO(Ll) 

HXSI-(  l.-XAD)*H0(L) +IAD*HO(H) 

HXXSI-(  1. -TAD)*HXO(L)  ♦TADmO(Ll) 

COMPUTE  THE  WATER  TEMPERATURE  AND  ITS  SPACE  DERIVATIVE. 

TXS I-Al*TO (  LI ) +A2*TO (  L)  ♦A3*TXO ( LI ) +A4*TX0 ( L) 

T(I)-((l.-0TAT2*ALPXSI)nxSI+DrAT2’»(ALP(I)*TA(KT)+ALPXSirrA(KTl))) 

*  /(l.+DXAT2*ALP(I)) 

TXXS  l«  B1  *TO  (  LI ) +82^0  (  L)  *83  *TXO  (  LI ) +B4*TX0  (  L) 

TX(I )-(  ( 1.  -0rAT2*(  VXXSI+ALPXSl)  )*TXXSI+OTAT2*(  (TA(KT1 )  -TXSI ) 

#  *ALPXSI*HXXSI/HXSI+(TA(KT)-T(1))*ALP(I)*HX(I)/H(I))) 

f  /(l.+0rAT2*(VX(I)+ALF(l))) 

CO  TO  200 

145  IF(CO(l)  .LE.  0.0  .AND.  TA(KT)  .GT.  0.0)  GO  TO  143 

LINEAR  INTERPOLATION  FOR  THE  FOLLOWING  QUANTITIES. 

VXXSI-(  1.  -TAU)*VXO(L)  +TAD«TX0(L1 ) 

HXS  I-(  1 . -TAU  )  *HO(  L) +IAU*HO(  LI ) 

HXXSX  l.-TAD)*HXO{L)+XAD’*HXO(Ll) 

COMPUTE  THE  ICE  CONTERATION  AND  ITS  SPACE  DERIVATIVE. 

CXS I-Al *CO ( U )  ♦A2*C0 (  L)  ♦A3*CXO (Ll ) +A4*CX0 ( L) 

C  ( 1 ) -CXS I-EPST2*  (TA(ICr  ) /H(  I ) +TA(rTl )  /  HXSI ) 

CXXSI-81*C0  (LI )  +B2*C0(L)  ♦83*CX0(L1 )  ♦B4*CX0(L) 

CX(I)-(< l.-DTAT2*VXXSI)*CXXSI 

#  +EPST2*(TA(KT)*HX{I)/H<I)**2+TA(KT1)*HXXSI/HXSI**2)) 

*  /(l.+DrAT2*VX(I)) 

GO  TO  200 

150  CONTINUE 

TRAJECTORY  OF  THE  CHARACTERISTIC  LINE  FALLS  ON  THE  D/S  BOUNDARY. 
INTERPOLATION  (»  TIME  IS  NEEKD. 

CQ-( VO( 1 ) -V( 1 ) ) /DELTAT 
CP»V(1)+V(I) 

CR— 2.*X(I) 

IF(CQ  .NE.  0.0)  GO  TO  155 
ETAS—CR/CP 
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GO  TO  156 

155  ETAS-(  -CP+(  CP*CP-4.  *CQ*CR  )**0. 5  )  /  (  2.  *CQ  ) 

156  XLAM-ETAS/ DELTA! 

IF(XLAM  .LE.  1.0)  GO  TO  160 
WRITE(6.400) 

400  FORMAT!/// ,5X,20C FAIL  TO  LOCATE  THE  STARTING  POINT') 

STOP  300 

160  ETAS2-ETAS*0.5 

IF(TO(I)  .LE.  0.0)  GO  TO  165 
C 

C  LINEAR  INTERPOLATION  FOR  THE  FOLLOWING  QUANTITIES. 

C 

162  TETA-(l.-XLAM)*T(l)+XLAM*TO(l) 

TAETA-(  1.-XLAM)*TA(1CT)+XLAM*TA(KT1) 

TTETA-( l.-XLAM)*TT(KT)+XLAM*TT(KTl) 

VETA-( 1. -XLAM)*V( 1)  +XLAM*VO( 1) 

VXETA“( 1. -XLAM)*VX(  1) +XLAM*7X0( 1) 

HETA-{  l.-XLAM)*H(  l)+XLAK*HO(l) 

HXETA-(  1.-XLAM)*HX(1)+XLAM*HX0<1) 

ALPETA-( l.-rLAM)*ALP(l)+XLAM*ALPO(l) 

C 

C  COMPUTE  THE  WATER  TEMPERATURE  AND  ITS  SPACE  lERIVATIVE. 

C 

T(I )-( <  1.  -ETAS2*ALPETA)*TETA+ETAS2*(ALPETA*TAETA+ALP(  1  )*TA(KT ) )  ) 

#  /(l.+ETAS2*ALP(I)) 

TXETA-( ALPETA*(TAETA-TETA )-TTETA ) / 7ETA 

TX(  I )  -(  (  1.  -ETAS2*  (  VXETA+ALPETA)  )*TXEIA+ETAS2*(  (TAETA-TETA) 

#  *ALPETA*HXETA/HETA+(TA(rr)-T(I )  )*ALP(I)»HX(  I )  /H(  I )  )  ) 

#  /(l.+ETAS2*(VI(I)+ALP(I))) 

GO  TO  200 

165  IF(CO(I)  .LE.  0.0  .AND.  TA(KT)  .GT.  0.0)  GO  TO  162 
C 

C  LINEAR  INTERPOLATION  FOR  THE  FOLLOWING  QUANTITIES. 

C 

CETA-( 1 . -XLAM)*C ( 1 ) 4XLAM*CO ( 1) 

CTETA-( l.-XLAM)*CT(Cr)+XLAK*CT{ETl) 

TAETA-( l.-XLAM)*TA(KT)+XLAM*TA(KTl) 

VETA»( l.-XLAM)*V( 1)  +X1AK*V0(1) 

VXETA»< l.-XLAM)*VX(l)+XLAM*VXO(l) 

HETA-( 1. -XLAM)*H( 1) +XLAM*HO( 1) 

HXETA=( 1.-XLAM)*HX(1)+XLAM*HX0(1) 

C 

C  COMPUTE  THE  ICE  CONCENTRATION  AND  ITS  SPACE  lERIVATIVE. 

C 

C  ( I )  =*CETA-EPS*ETAS  2*(TA  (XT  ) /H(  I ) +TAETA/HETA) 

CXETA—  (  EPS*TAETA  /HETA+CTETA)  /  VETA 
CX(I)-( l.-ETAS2*VXETA)*CXETA 

#  ♦EPS*ETAS2*(TA(KT)*HX(I)/H(I)**2+TAETA*HXETA/HETA**2) 

200  CONTINUE 

DO  210  I-1,N 
HO(I)-H(I) 

HXO(I)-HX(I) 

V0(I)-V(I) 

VXO(I)-VX(I) 

T0(I)-T(I) 

TXO(I)-TX(I) 

CO(I)-C(I) 

CX0(I)-CX(I) 

ALPO(I)-ALP(I) 

210  CONTINUE 
RETURN 
500  FORMAT! 

#'  <  THERMAL  CONDITION  >',// 

#'  Energy  exchange  coeff.  at  the  air-water  interface  - 

#  F6.2,'  cal/cm**2.day.C',/ 

#'  Density  of  water  “  ',F5.2,'  g/ciB**3',/ 

#'  Specific  heat  of  water  “  ',F5.2,'  cal/g.C',/ 

#'  Latent  heat  of  fusion  of  ice  “  ',F6.2,'  cal/g',//) 

END 
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SDB ROUTINE  DIFEX(F,  PX.H) 

COMFOTE  THE  SPACE  BERIVATIVL  USING  THE  WEIGHTED  DIFFERENCE. 

COMMON  /C13/  DELX(34),I(34) 

DIMEN3I0..  F(1>,FX(1} 

Nl-N-l 

FX( 1) -( F( 2)-F( 1} )/DELX(  D 
FX(!i)-(F(N)-F{Nl))/DKLX(Nl) 

CO  10  I«2,H1 
Il-I-l 

51- (F(I)-F(I1))/DEI.X(I1) 

52- ( F( I+l ) -F(I ) ) /DELX( I ) 

10  rX(I)-(Sl*DELX(I)+S2*DELX(Il))/(DELX(H)+IlELX(I)) 

RETURN 

END 


SUBROUTINE  DIFFTCF,  FT.NT) 

COMPUTE  THE  TIME  lERIVATIVE  USING  THE  CENTRAL  DIFFERENCE. 

COMMON  /C09/  GRAV.SWI.DELTAT,ALFHA,BETA,THETA 
DIMENSION  FCJ.FTd) 

NTl-NT-i 

FT(1)-(F(2)-F(1))/DELTAT 
IT ( NT ) -< F( NT) -F ( NTl ) ) / DELTAT 
DO  10  J-2,NT1 

10  FT(J)-(F(J+1)-F(J-1))/(2.*DELTAT) 

RETURN 

END 


jUBROnTINE  AIRTMP 

DEFINE  THE  BOUNDARY  CONDITION  FOR  TEE  AIR  TEMPERATURE. 

COMMON  /C06/  N.NT.TIME.I.KCiMfTl.LE 
COMMON  /C09/  GBAV,SWI,mTAT, ALPHA,  BETA, THETA 
COMMON  /C16/  TA(200),TTDB(200),TO<34),CrUB(200),CO(34) 
READd,*)  CONTA 
DO  10  .<->l,HNTl 
10  TA(J) -CONTA 
RETURN 
END 


SUBRCRITINE  TICOND 

DEFINE  THE  INITIAL  CONDITION  FOR  THE  WATER  TEMPERATURE. 
COMMON  /C06/  N,NT,TIME,I,Kr,MNTl,LE 

COMMON  /C16/  TA(200),TTUB(200),T0(34),CrUB(200),CO(34) 
READd,*)  CONTO  ' 

DO  10  I-1,N 
10  TO(I)-COHTO 
RETURN 
END 
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SUBROUTINE  TUPBC 
C 

C  DEFINE  THE  U/S  B.C.  FOR  THE  WATER  TEMPERATURE. 

C 

COMMON  /C06/  N.NT.TIME.I.KT.MHTI.LE 

COMMON  /C09/  CRAV, SHI, DELTAT, ALPHA. BETA, THETA 

COMW5N  /C16/  TA(200),TTUB(200),TO(34),CTUB(200J,CO(34) 

READd,*)  CONTUB 

DO  10  J»1,MNT1 

C  10  TTUB(J)-5.-5,*SIN(2.*3.141592*(J-1)*DELTAT/36000.) 

10  TTUBCJ)  “CONTUB 
RETURN 
END 

SUBROUTINE  CICOND 
C 

C  DEFINE  THE  INITIAL  CONDITION  FOR  THE  ICE  CONCENTRATION. 
C 

COMMON  /C06/  N.NT.TIME.I.KI.MNTI.LE 
COMMON  /C16/  TA(200),TTCB{200),TO(34),CrUB{200>.CO(34) 
DO  10  1-1 ,N 
10  C0(l)-0.0 
RETURN 
END 


SUBROUTINE  CUPBC 
C 

C  MFINE  THE  U/S  B.C.  FOR  THE  ICE  CONCENTRATION. 

C 

COMMON  /C06/  N.NT.TIME.I.KT.MNTl.LE 
COMMON  /C09/  GRAV.SWI.IELTAT.ALPHA.BETA.THBTA 
COMMON  /C16/  TA(200),TTUB(200),TO(34),CrDB(200}.CO(34) 
DO  10  J-1,KNT1 
10  CTDB(J)-0.0 
RETURN 
END 


SUBROUTINE  ICECVR 
C 

C  SIMULATION  OF  THE  ICE  COVER  FORMATION 
C 

COMMON  /COl/  A(34).B<34).BW(34),H(34>,R(34),S(34) , 
#  CK(34),m(34).CKST(34),Z(34) 

COMMON  /C02/  Q(34),V(34),T(34) 

COMMON  /C06/  N.NT.TIME.I.KT.MNTl.LE 

COMMON  /C09/  GRAY, SHI, DELTAT, ALPHA, BETA, THETA 

COMMON  /C13/  DELX(34).X<34) 

COMMON  /C15/  ICEFLG(34),T(34).C(34),TH(34) 

1F(KT  .GT.  1)  GO  TO  5 
Nl-N-1 

READd,*)  FRMAX,FRMIN,TBFLOE,ALPS,EI,BP,SWI 
EC-EP+(1.-EP)*EI 

WRITE(6,100)  FRMAX.FRMIN.TEFLOE.ALPS.EI.EP.EC.SHI 
CONA-(  2.  *(  l.-EC)*d.-SHI)  )**0.5 
5  CVOL-ALPS*C(LE)*q(LE)*M:LTAT 
IF(CVOL  .LE.  0.0)  RETURN  ' 

FR- V(  LE)  /  (  GRAV*H(  LE)  )  **0. 5 
IP(FR  .  GT.  FRMAX)  GO  TO  60 
IF(FR  .GT.  FRMIN)  GO  TO  10 
C 

C  FROUDE  NO.  IS  LESS  THAN  MINIMUM  FROUDE  NO. 
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THE  THICKNESS  OF  THE  ICE  COVER  HItL  BE  EQDAL  TO  THE  THICKNESS 
OF  ICE  ElOE. 

THMBX-THFIOE 
GO  TO  35 

FHOOIE  NO.  IS  BETWEEN  THE  MAXIMUM  AND  THE  MINIMUM. 

COMPUTE  THE  ICE  COVER  THICKNESS  BASED  ON  THE  ICE  JAM  THEORY. 

10  KER-l 
J^*0.0 

20  F-H*(W*W-1.0)+FR/CONA 
FD-3.0*W*V-1.0 
H-H-F/FD 

IF(ABS(F}  .LT.  l.B-6  .OR.  ITER  .GE.  101  GO  TO  30 
ITER-ITER+1 
GO  TO  20 

30  THMAX-W*W*H(LE) 

35  DTH-CV0I,/(DEUC(LE)*(B(LE)+B(LE+1))*0.5) 

IF(TH(l.E)+D(rH  .LT.  THMAI)  GO  TO  50 
DVOL-CVOL-CVOL*(THMAX-IH(  US ) ) /DTH 
TH(LE}-TEMAX 
C(LE-H)-0.0 

IF(LE  .EQ.  HI)  TH(N)-TH(LE) 

LE-LE-1 

IP(LE  .GT.  0)  GO  TO  40 
WRITE(6,36) 

36  F0RMAT( ////,'  ********  ICE  COVER  REACHES  THE  U/S  *********' j 

#  '  irirkirtrirk*  PROGRAM  TERMIHATED  *********',////) 

STOP 

40  TH(LE)-DV0L/(DELX(LE)*(B(LE)+B(LE+1)>*0.5) 

RETURN 

50  TH(LE)-TH(LE)+DrH 
C(LE+l)-0.0 

IF(LE  .BQ.  Nl)  TH(N)-TH(LE) 

RETURN 

60  IF(LE  .KE.  Nil  GO  TO  65 
WRITE(6,631  FR 

63  FORMATC/////,'  *******  NO.  OF  D/S  -  ',F5.2, 

#  'IS  GREATER  THAN  THE  IHMAX  :  NO  ICE  COVER  CAN  BE  FORMED 

#  *****',/////) 

STOP 

C 

C  FROUDE  NO.  IS  GREATER  THAN  MAXIMUN  FROUDE  NO. 

C  FIND  THE  NEAREST  D/S  SECTION  WHERE  THE  FROUDE  NO.  IS  LESS  THAN 

C  MAXIMUM  FROUDE  NO.  AND  THEN*  ACCUMULATE  THE  BRAZIL  ICE  AT  THAT 

C  SECTION. 

C 

65  LEl-LE+1 

DO  70  K-LEl  ,N1 
FR-V(K)/(GRAV*H(K))**0.5 
IFIFR  .L£.  FRMAX)  GO  TO  SO 
70  CONTmUE 

80  0IH-CVOL/{DELX(K)*(B(K)+B(K+U)*O.5) 

TH(K)-TH(K)+DTH 
IF(K  .BQ.  Nil  TH(H)-TH(K) 

RETURN 
100  FORMATC 

#'  <  ICE  CONDITION  >',// 

#'  Maxinun  Fr  No.  at  vhich  the  ice  cover  can  not 

# 

#'  MiniiBaiB  Fr  No.  for  the  fomation  of  ice  cover 
#'  The  thickness  of  ice  floe  »  ',F5.3,'  ■',/ 

#'  The  ratio  of  the  surface  ice  discharge  to  the 
#'  discharge  •  ',F4.2,/ 

#'  Porosity  of  individual  ice  floes  ••  ',F5.3,/ 

#'  Porosity  of  the  ice  accumulation  between  ice  floe*  “  ',F5.3, 


progress  •> 
F5.2,/ 

-  ',F5.2,/ 


total  ice' I 
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Porosity  of  the  ice  sccumulstion  -  ',P5.3,/ 
#'  Specific  gravity  of  the  ice  “  ',P5.3,//) 

END 


SDBRCHITINE  THICON 

COMMON  /C06/  N,NT,TIME,I.KT,MNT1.LE 
COMMON  /CIS/  ICEnx;<34).T(34),C(34),TH(34) 
DO  10  I-1,N 
10  TH(I)-0.0 
RETURN 
END 


SUBROUTINE  CONDUC 

COMMON  /C06/  N,HT,T1ME,I,KT.MNT1,LE 

COMMON  /C09/  GREV, SHI, DELT4T, ALPHA, BET&.TBETA 

COMMON  /CIS/  ICEELG(34),T(34),C(34),TH(34) 

COMMON  /C16/  TA(200),TTUB(200),TO(34),CTUB{200>,CO(34) 

COMMON  /C17/  ROU,BtfA,CLAIHT 

IF(KT  .GT.  1)  GO  TO  10 

READd,*)  HU.TCI 

HR1TE(6,200)  HU,TCI 

ROni-ROU*$HI 

HIA-HIA*10000.  /(  24.  *3600. ) 

TCI-TCI*100. 

ETRC-DELTAT/(BOUI*CLATHT) 

10  DO  100  I-LE,M 

IF(TH(I)  .LE.  0.0)  GO  TO  100 
TS-TAiKT )  *aiA*TH(  I)/ (  HU*TH(  I ) +TCI ) 

IF(TS  ,GT.  0.0)  GO  TO  SO 
C  DrrH«-(TS*TCI/TH(l)+HHI*T(I))*DTRC 

DTH=-TS *TCI /TH ( I ) *DTRC 
TH(I)-TH(I)+DTH 
GO  TO  100 

C  DTH— (HU*TA(KT)+HHI*T(I))*DTRC 

SO  DTH— HU*TA(KT)*0TRC 
TH(I)»TH(I)+OrH 
IF(TH(I)  .LT.  0.0)  TH(l)-0.0 
100  CONTINUE 
RETURN 
200  FORMAT ( 

f'  Energy  exchange  coeff.  at  the  air-ice  interface  •  ' 
#  F6.2,'  cal/ciB**2.day.C',/ 

#'  Thermal  conductivity  of  ice  »  ',F7.4,''  cal/cm.sec.C 
END 
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cially  large  ones,  the  most  appropriate  method  involves  controlled  ice-cover  formation  through  the  regulation  of  river  flow. 
Flow  discharge  and  stages  would  be  regulated  by  controlling  the  flow  releases  from  reservoirs  and  flow  stages  at  river 
dams,  such  that  optimal  flow  conditions  prevail  for  rapid  formation,  and  subsequent  maintenance,  of  an  accumulation  ice 
cover  over  river  reaches  in  which  potentially  large  amounts  of  frazil  ice  may  grow.  Accumulation  covers  would  be  formed 
of  frazil  ice  pans  and  floes  and,  if  appropriate,  broken  ice  conveyed  from  upstream.  Existing  dams,  augmented  where  need¬ 
ed  by  navigable  ice  booms,  could  serve  as  retention  structures  for  the  development  of  accumulation  covers.  A  preliminary 
indication  of  the  feasibility  of  this  method  for  controlling  ice-cover  formation  on  stage-regulated  pools  of  the  Ohio  River  is 
assessed  through  the  use  of  a  numerical  model  that  simulates  ice-cover  formation  from  frazil  ice.  It  is  found  that  this  ap¬ 
proach  holds  promise  for  mitigating  jam  occurrence,  although  its  implementation  necessarily  entails  management  of  flow 
through  major  portions  of  the  Ohio  River.  The  results  of  the  study  are,  to  a  limited  extent,  generalized  to  other  rivers. 
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